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Abstract 


This  report  describes  a  package  which  simulates  the  dynamic  radar  signatures  a  radar 
receiver  would  expect  when  tracking  a  moving  target  such  as  an  aircraft  or  a  missile 
The  package  consists  of  several  MATLAB  codes.  The  basic  engagement  simulated  is 
a  radar,  such  as  seen  on  a  fully  active  missile  seeker,  engaging  a  target  such  as  an 
aircraft.  The  static  scattering  characteristics  of  the  target  are  required  as  an  input 
to  the  package.  The  package  then  simulates  the  radar-target  engagement  dynamics 
modeling  the  target  motion  and  the  physical  target  vibrations  which  effect  the  signa¬ 
ture  the  radar  receives.  The  outputs  of  the  package  are  the  scattered  electric  field, 
the  Doppler  shift,  and  the  angular  tracking  error,  or  glint,  all  as  seen  by  the  radar. 
This  package  is  meant  to  aid  a  hardware-in-the-loop  simulation  of  a  radar-target 
engagement. 
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1.  Introduction 


This  report  describes  a  package  which  simulates  the  dynamic  radar  signatures  a 
radar  receiver  would  expect  when  tracking  a  moving  target  such  as  an  aircraft  or  a 
missile.  By  dynamic,  it  is  meant  that  the  radar  and/or  target  are  physically  moving. 
The  package  consists  of  several  MATLAB  codes. 

This  work  was  performed  at  Australia’s  Defence  Science  and  Technology  Organi¬ 
sation  (DSTO),  in  Salisbury,  South  Australia,  Australia,  under  the  United  States’  Air 
Force  Office  of  Scientific  Research  Scientist  and  Engineer  Exchange  Program,  from 
February  1,  1996  to  May  10,  1996,  by  the  author.  The  office  within  DSTO  where  the 
work  was  performed  is  The  Maritime  and  Aeronautical  Research  Laboratory,  Guided 
Weapons  Division,  Radio  Frequency  Seekers  Group.  The  supervisor  at  DSTO  was  Dr. 
Noel  Martin,  Principal  Research  Scientist  and  Head  of  the  Radio  Frequency  Seekers 
Group. 

The  basic  engagement  simulated  is  a  mono-static  radar,  such  as  seen  on  a  fully 
active  missile  seeker,  engaging  a  target  such  as  an  aircraft.  The  radar  is  assumed 
to  be  mono-static  (co-located  receive  and  transmit  antenna).  The  static  (no-motion) 
scattering  characteristics  of  the  target  are  required  as  an  input  to  the  package.  These 
are  available  for  a  representative  target,  the  C-29  transport  aircraft.  The  package 
uses  a  simplified  model  for  the  static  scattering  from  the  target,  the  2-D  Damped 
Exponential  Model.  This  model  is  sometimes  called  the  multi-point  scattering  model 
and  is  described  is  several  references  [1,  2,  3]. 

The  package  then  simulates  the  radar-target  engagement  dynamics  modeling  the 
target  motion  and  the  physical  target  vibrations  which  will  effect  the  signatures  the 
missile  receives.  The  outputs  of  the  package  are  the  scattered  electric  field,  the 
Doppler  shift,  and  the  angular  tracking  error,  or  glint,  all  as  seen  by  the  radar. 
This  package  is  meant  to  aid  a  hardware-in-the-loop  simulation  of  a  radar-target 
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engagement.  Figure  1.1  shows  a  flow  diagram  outlining  the  development  of  this 
simulation  package. 

The  notational  conventions  followed  in  this  report  are  that  vectors,  or  quantities 
with  a  magnitude  and  direction,  are  so  noted  by  the  “vec”  notation,  for  example, 
the  vector  Rt.  Unit  vectors,  or  vectors  with  unit  magnitude,  as  so  designated  by 
a  “hat”,  for  example  the  unit  vector  Dt.  Also,  complex  quantities  are  denoted  by 
boldtype,  for  example,  the  complex  number  Px-  The  imaginary  basis  is  denoted  by 
j,  i.e.  j  =  Scalar  quantities  do  not  have  the  “vec”  sign  and  are  not  boldtyped. 

Throughout  this  report,  the  term  radar  is  used  to  describe  what  is  tracking  the 
target.  The  radar  can  be  a  conventional  radar  as  well  as  a  missile  seeker.  Also,  the 
radar  can  be  stationary  or  on  a  moving  platform  such  as  an  aircraft.  The  term  target 
describes  what  is  being  tracked  by  the  radar.  The  target  can  be  an  aircraft,  or  any 
flying  or  stationary  target.  In  this  report,  the  radar  is  always  called  the  radar  while 
the  target  is  usually  called  the  target,  but  is  also  called  the  aircraft  in  many  places, 
since  the  target  shown  in  the  simulations  is  a  C-29  transport  aircraft. 

This  report  is  organized  as  follows.  Section  2  covers  the  model  used  for  the  target 
scattering  and  defines  the  dynamic  scenario  between  radar  and  target.  Section  3 
defines  the  output  quantities  of  the  simulations.  Section  4  shows  the  results  of  a 
typical  simulation.  Section  5  describes  the  MATLAB  codes  used  to  implement  the 
simulations.  Finally,  Section  6  concludes  the  report  and  provides  recommendations 
for  further  research  in  the  modeling  of  dynamic  radar  data. 
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C-29  MODEL 

2-D  Damped  Exponentials 

STATIC  DATA 


FLY”  C-29  VS.RADAR 
Via  Simulation 
profile,  m 
ADD  MOTION 


scattered.field.  m 
Quantities  Radar  Receiver 
would  see  on  typical  C-29 
radar  encounter 


Figure  1.1:  Simulation  development  methodology.  First,  simplified  static  scattering 
models  are  developed.  Then  motion  compensation  is  added  to  the  models.  Finally, 
the  quantities  of  interest  are  calculated. 
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2.  Radar- Target  Engagement  Geometry,  Dynamics,  and  Target 

Scattering 

This  section  describes  the  radar-target  dynamics  and  the  scattering  model  used 
for  the  target.  Along  with  the  descriptions  of  the  dynamic  situations  and  target 
scattering  calculations,  examples  are  used  to  demonstrate  the  theory. 

2.1  Basic  Flight  Path 

The  basic  engagement  geometry  is  shown  in  Figure  2.1,  The  variables  used  in  the 
engagement  geometry  are  described  in  Table  2.1.  Please  note  that  all  units  are  in  the 
metric  system,  i.e.,  all  angles  are  in  degrees,  all  frequencies  are  in  hertz,  all  distances 
in  meters,  all  velocities  in  meters  per  second,  all  powers  in  watts,  and  all  gains  are  in 
numeric  (not  in  dB). 

For  the  simulations  it  is  possible  that  both  the  radar  and  target  are  moving.  In  a 
real  situation,  of  course,  this  is  the  case,  but  for  the  simulations  performed  here  it  is 
assumed  that  the  radar  is  fixed  at  the  origin  of  the  coordinate  system.  Thus,  a  radar 
fixed  coordinate  system  is  assumed  and  Ro  =  (0, 0, 0).  Also,  it  is  assumed  that  the 
radar  illuminating  the  target  is  continuous- wave  (CW),  having  a  100%  duty  cycle. 
If  a  pulsed  radar  is  desired,  the  results  of  the  simulations  can  easily  be  converted  to 
account  for  this  [4] .  The  angles  AZ  and  EL  are  relative  to  the  waterline,  or  centerline, 
of  the  target,  defined  by  the  direction  vector  {Pt). 

The  simulations  are  set  up  so  that  at  a  given  instant  in  time,  t,  the  target  is 
at  a  given  position,  Rt,  traveling  in  the  Dt  direction  at  velocity  Vt.  The  target  is 
assumed  to  be  pointing  in  the  Pt  direction.  By  pointing,  it  is  meant  that  the  physical 
centerline,  or  waterline,  of  the  target  is  in  the  Pt  direction.  For  now,  we  assume  the 
Dt  =  Pt.  In  Section  2.3,  roll,  pitch,  and  yaw  angles  are  introduced  to  allow  Dt  and 
Pt  to  take  on  different  values.  Next,  all  of  the  simulation  outputs  are  calculated. 
Then  time  is  incremented  by  dt,  the  target  (aircraft)  moves  the  appropriate  distance. 
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Figure  2.1:  Radar  Target  Geometry. 

and  the  outputs  are  calculated  again,  and  so  on.  Also,  the  normalized  time  variable, 
n,  increments  as  n  =  1, 2, . . .  where  there  are  ti  time  steps  in  a  simulation.  Thus, 
the  actual  time  t,  in  seconds,  is  given  by  t  =  (n  -  l)dt. 

The  MATLAB  code  which  runs  the  entire  simulation  and  calls  all  of  the  MATLAB 
functions  is  called  profile.m.  This  is  a  MATLAB  script  file  and  is  commented  to 
show  where  the  various  functions  are  performed. 

In  the  simplest  case,  the  target  is  fiying  straight  and  level  with  no  perturbations 
(vibrations)  in  position  or  orientation.  In  reality  this  is  not  true.  The  target  has 
perturbations  in  both  position  (translational)  and  orientation  (rotational)  from  the 
intended  flight  path.  Although  these  perturbations  are  very  small  relative  to  the 
size  of  the  target,  they  are  the  driving  factor  behind  angle  error  (glint)  and  Doppler 
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Table  2.1:  Simulation  Variables. 


Variable 

Description 

Units 

/ 

radar  frequency 

hertz 

Ptr 

transmitter  power 

watts 

Gtr 

transmitter  antenna  gain 

numeric 

Vt 

target  velocity 

meters  per  second 

Di 

target  direction  vector  (x,y,z) 

normalized 

Pi 

target  pointing  vector  (x,y,z) 

normalized 

Rt 

target  position  vector  (x,y,z) 

meters 

Ro 

radar  position  vector  (x,y,z) 

meters 

dt 

time  step  for  simulations 

seconds 

ti 

number  of  iterations 

integer 

t 

time 

seconds 

n 

normalized  time 

integer 

AZ 

azimuth  angle 

degrees 

EL 

elevation  angle 

degrees 

shift.  Thus,  these  perturbations,  or  vibrations,  are  modeled  here  and  described  in 
Section  2.5.  Note,  as  stated  before,  the  target  can  be  any  flying  object,  but  it  may 
also  be  referred  to  as  the  aircraft  since  the  aerodynamic  models  assume  an  aircraft. 

The  scattering  models  developed  in  the  next  section  apply  to  both  the  far-fleld 
and  the  near-field  scattering  from  a  target.  Thus,  the  simulations  developed  can  be 
applied  to  large  radar-target  range  —  R6)  scenarios  (far-fleld)  or  for  small  radar- 
target  range  (near-field  or  end  game)  scenarios.  Also,  in  the  course  of  a  simulation, 
both  far-fleld  and  near-field  situations  can  exist,  and  the  simulations  can  accurately 
predict  the  scattered  field,  Doppler  shift,  and  angular  error  throughout  this  scenario. 

Next,  in  Section  2.2,  the  aircraft  scattering  model  will  be  described.  This  simple 
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mathematical  model  calculates  the  scattered  field  from  any  radar  target,  such  as  an 
aircraft,  as  a  function  of  polarization,  frequency,  and  aspect  angle  for  a  static  situation 
(stationary  radar  and  target). 

2.2  Aircraft  Scattering  Model 

In  order  to  calculate  the  scattered  field  the  radar  receives,  the  mono-static  radar 
scattering  characteristics  of  the  aircraft  must  be  modeled.  This  section  describes 
how  the  static  (no-motion)  scattering  response  of  a  target  is  calculated.  The  model 
used  relies  on  the  scattering  center  approach  to  scattering,  which  assumes  that  the 
high  frequency  scattering  from  an  electrically  large  target  such  as  an  aircraft  can  be 
approximated  by  a  sum  of  scattering  from  a  finite  number  of  predominant  scattering 
centers  [5,  6,  7,  2].  The  single-polarization  scattering  from  a  complicated  target  is 
thus  approximated  as 

=  (2.1) 

where  S{f,0)  is  the  scattering  coefiicient  for  the  entire  target  and  can  be  further 
defined  as 

S(/,^)  =  ^Es^(/,0)j  ,  (2.2) 

where  {f,e)  is  the  scattering  coefficient  for  the  q-th  scattering  center,  {f,e)  is 
the  scattered  electric  field  (frequency  domain)  as  seen  at  the  radar  receiver,  (/,  9) 
is  the  electric  field  transmitted  from  the  radar,  and  r  is  the  distance  from  the  radar 
to  the  aircraft  (target)  phase  center.  The  scattering  coefficient,  and  the  scattered 
fields  are  all  complex  numbers  and  also  functions  of  angle  and  frequency.  They  are 
considered  frequency  domain  representations  of  the  electric  fields. 

The  term  aircraft  phase  center  refers  to  a  point  on  or  near  the  physical  target 
center  to  which  all  of  the  phases  of  (/,  0)  are  referenced.  Then  r  is  the  distance 
from  this  point  to  the  radar  receive  antenna.  The  angle  ^  is  a  generic  angle  which 
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defines  radar-target  orientation  as  referenced  from  the  target.  However,  for  all  of 
the  simulations  in  this  report,  Q  is  the  azimuth  angle  (AZ)  relative  to  the  aircraft, 
as  shown  in  Figure  2.2.  The  frequency  /  is  the  frequency  of  the  CW  electric  field 
transmitted  by  the  radar. 

Also,  this  is  considered  a  single  polarization  version  of  scattering,  meaning  that  the 
polarizations  of  the  electric  fields  are  user  definable,  but  fixed  for  a  specific  model. 
In  all  of  the  simulations,  the  hh  polarization  (horizontally  polarized  transmit  and 
receive  antenna  on  the  radar)  is  chosen,  and  the  C-29  transport  aircraft  is  modeled. 
If  different  transmit  or  receive  polarizations  are  required,  then  a  different  scattering 
coeflicient,  S[f,d),  will  result.  In  the  bi-static  case  (different  transmit  and  receive 
antenna  locations),  then  a  different  scattering  coefficient  must  be  used. 

For  all  the  simulations  performed  here,  the  phase  of  the  transmitted  field  is  as¬ 
sumed  to  be  zero  at  the  transmit  antenna,  while  the  magnitude  is  determined  by  the 
radar  power,  Ptr,  and  the  radar  antenna  gain  Gtr.  It  is  assumed  in  all  simulations 
that  the  radar  antenna  is  bore-sighted  directly  at  the  target.  Further,  it  is  assumed 
that  the  transmission  media  between  the  target  and  the  radar  is  lossless,  isotropic 
free  space.  Thus,  the  propagation  losses  between  the  radar  and  the  target  (two-way) 
are  taken  into  account  by  the  term. 

The  r  scattering  coefficients,  s.^  if,0),  contain  all  the  target  information.  In  this 
work,  the  scattering  coefficients,  which  are  functions  of  frequency  (/)  and  radar- 
target  aspect  angle  {$),  are  modeled  with  a  2-D  Damped  Exponential  Model.  The 
2-D  Damped  Exponential  Model  is  given  by 

r 

d(m,  n)  =  J2  a7Px,™Py/  (2.3) 

7=1 

where 

Px.,  =  a;-pole,  ai-component  of  2-D  exponential  (complex  number) 

Py7  =  Tfh  y-pole,  y-component  of  2-D  exponential  (complex  number) 
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=  7th  amplitude  coefficient  (complex  number) 
r  =  number  of  scattering  centers  (integer) 
m  =  0, 1, ...  M  —  1;  (normalized  frequency) 
n  =  0, 1, ...  iV  —  1;  (normalized  angle). 

This  is  a  2-D  extension  of  Prony’s  Model  [2].  This  model  assumes  the  radar  target  is 
comprised  of  P  scattering  centers.  The  7-th  scattering  center’s  scattering  behavior, 
for  a  single  transmit  and  receive  polarization,  is  modeled  by  the  x-pole,  y-pole,  and 
amplitude  coefficient  triple,  {px^, Py^, a^}. 

It  is  emphasized  that  this  is  a  2-D  frequency  domain  model.  A  2-D  Inverse  Fourier 
Transform  of  data  in  this  domain  yields  data  in  the  ISAR.  Image  domain  (or  range 
domain).  See  Appendix  A  for  a  detailed  description  of  this  relationship  and  how  an 
ISAR  image  is  formed.  In  this  report,  only  ISAR  image  plots  are  shown,  not  plots  of 
data  in  the  2-D  frequency  domain,  as  plots  of  data  in  the  frequency  domain  contain 
very  little  useful  visual  information. 

The  magnitude  of  the  x-component  damped  exponential,  |px^|,  determines  the 
dispersion  of  the  7-th  scattering  center  in  the  ^-coordinate  of  the  transform  domain 
(down-range  in  the  ISAR  image).  By  dispersion,  it  is  meant  how  wide  the  response  of 
the  7-th  scattering  center  is  in  the  down-range  direction  of  the  ISAR  image  (similar 
to  a  point  spread  function  [8]).  The  angle  of  the  x-component  damped  exponential, 
determines  the  location  of  the  7-th  scattering  center  in  the  down-range  direc¬ 
tion.  Similarly,  the  magnitude  of  the  ^-component  damped  exponential,  |py^|,  deter¬ 
mines  the  dispersion  of  the  7-th  scattering  center  in  the  y-coordinate  of  the  transform 
domain(cross-range  in  the  ISAR  image)  while  the  angle  of  the  y-component  damped 
exponential,  /Py.,,  determines  the  location  of  the  7-th  scattering  center  in  the  cross¬ 
range  direction  of  the  ISAR  image.  For  this  work,  the  x-poles  correspond  to  the 
down-range  dimension  and  the  y-poles  correspond  to  the  cross-range  dimension  in  an 
ISAR  image  (see  Figure  2.5). 
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Relating  the  2-D  Damped  Exponential  Model,  Equation  (2.3),  to  the  model  in 
Equation  (2.1),  we  can  see  that 


S{f,e)  =  d(m,n). 


(2.4) 


or 


and  that 


r 


r 


E^(/.«)  =  Ea7P.,“Py,". 

7=1  7=1 


s,  (/,«)=  a, p»,”Py,* 


(2.5) 


(2.6) 


for  7  =  1,2,  ...,r.  It  must  be  noted  that  the  normalized  frequency  m  is  directly 
related  to  the  actual  frequency  /,  and  the  normalized  angle  n,  is  directly  related  to 
the  aspect  angle  to  the  aircraft  6  (see  Equation  (2.11)).  The  angle  6  is  directly  related 
to  the  azimuth  and  elevation  angles  at  which  the  aircraft  is  being  illuminated  from. 
For  all  of  the  C-29  data, 

e  =  Az. 


Thus,  all  of  the  C-29  data  are  azimuth  cuts  consistent  with  the  definition  of  AZ  shown 
in  Figure  2.2.  Note  that  the  definition  of  AZ  in  the  figure  is  negative  azimuth  (-AZ). 
Positive  AZ  values  go  in  the  clockwise  direction  looking  from  above  the  aircraft.  Also 
shown  in  Figure  2.2  is  the  target  direction  vector,  Pt,  as  defined  in  Table  2.1. 

From  the  geometry  of  the  situation,  the  values  of  AZ  and  EL  from  the  aircraft  to 
the  radar  can  be  determined.  From  these  angles,  values  of  m  and  n  can  be  found. 
The  values  of  m  and  n  depend  upon  how  the  values  of  {Px.,,  Py.^, 87}  were  originally 
determined. 

Consider  samples  of  the  scattering  coefficient  S{f,9)  taken  at  frequencies  /  = 
/O)  fly-,  fu-i  and  angles  0  —  60,  Oi,..  .,6^-1.  It  is  assumed  that  the  frequencies 
and  angles  are  evenly  spaced  at  intervals  of  df  and  da,  respectively.  If  the  data  is 
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Figure  2.2:  Definition  of  azimuth  (AZ)  and  elevation  (EL)  angles  for  the  target.  The 
AZ  angle  shown^  actually  negative  AZ  (-AZ).  Note  that  for  all  the  C-29  data,  0=k'L. 
Also  note  that  Pt  is  the  direction  the  aircraft  is  pointing  in. 

formed  into  a  matrix,  an  M  x  iV  data  matrix  called  D  results,  and  is  given  in  by 

S(/o,^o)  S(/o,0i)  •••  S(/o,%_i) 

S(/i,0o)  S(/i,0i)  ...  S(/i,^^_x) 

S{fM-X,eo)  S(/m-1,^i)  ••• 

The  total  frequency  and  angular  bandwidths  of  the  data  are  defined  by 

AF  =  fM-i-fo  (2.8) 

A0  =  6n-i  —  ^0  (2.9) 

As  stated  before,  from  this  data,  an  Inverse  Synthetic  Aperture  Radar  (ISAR)  Image 
can  be  formed  using  the  2-D  Inverse  Fourier  Transform.  For  a  brief  description  of 
radar  imaging,  please  see  Appendix  A  or  consult  [8,  9]. 

From  this  raw  radar  data  contained  in  the  matrix  D,  the  2-D  Damped  Exponen¬ 
tial  Parameters,  {Px.,,  Py,,  are  estimated  using  the  2-D  TLS-Prony  Technique. 

This  technique,  described  in  [3],  takes  the  raw  radar  data  and  finds  the  “best”  pa¬ 
rameters  in  a  Total  Least  Squares  sense.  Once  these  parameters  are  found,  they  can 
be  used  over  and  over  again  to  predict  the  scattering  from  a  target  for  a  range  of 
frequencies  and  angles.  The  raw  ISAR  data  in  D  is  considered  “unfocused”  and  is 
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used  directly,  without  any  focusing  (see  Appendix  B),  to  estimate  the  2-D  Damped 
Exponential  parameters. 

It  must  be  noted  that  the  model,  described  by  is  only  valid  over 

a  range  of  m  and  n  related  to  the  2-D  angular  and  frequency  range  of  the  original 
raw  ISAR  data,  AF  and  A0,  from  which  the  model  was  developed.  Using  the  model 
outside  of  these  limits  is  warranted,  but  only  for  a  limited  number  of  frequencies 
and  angles.  A  good  rule  of  thumb  is  to  use  the  2-D  Damped  Exponential  Model 
at  frequencies  between  (/o  -  AF)  and  (fu-i  +  AF)  and  angles  between  (Oq  -  A0) 
and  -f  A0).  This  corresponds  to  values  of  normalized  frequency  m  between 

(— M  -I- 1)  and  (2M  —  1)  and  values  of  normalized  angle  n  between  (— y)  and 
For  more  information  on  this  topic,  which  is  called  data  extrapolation,  see  [10].  Also, 
the  normalized  frequency  and  angular  variables,  m  and  n,  are  related  to  the  actual 
variables,  /  and  6,  as 


m 

n 


f-fo 


df 


o-ep 

da 


+  iV/2. 


(2.10) 


Now  that  the  modeling  procedure  has  been  described,  an  example  of  how  it  was 
applied  to  some  actual  static  radar  data  is  shown.  For  a  more  detailed  explanation 
of  the  2-D  Damped  Exponential  Model  and  its  application  to  radar  modeling,  please 
consult  [1,  2,  3,  11]. 

2.3  Damped  Exponential  Modeling  of  C-29  Scale  Model 


The  MATLAB  code  which  calculates  the  2-D  Damped  Exponential  parameters 
for  the  C-29  data  is  called  RUN.m.  The  MATLAB  code  which  contains  the  model 
parameters  for  the  C-29  aircraft  is  called  C29model.m.  The  raw  C-29  data  was 
obtained  from  a  compact  Inverse  Synthetic  Aperture  Radar  (ISAR)  range,  and  the 
parameters  of  the  raw  ISAR  data  are  shown  in  Table  2.2.  The  C-29  aircraft  measured 
was  a  high  quality  one-third  scale  model  of  the  actual  C-29  aircraft.  An  outline  of 
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the  scale  model  is  given  in  Figure  2.3.  For  more  information  on  the  C-29  aircraft 
consult  jane’s  All  the  World’s  Aircraft  or  [12].  The  one-third  scale  C-29  model  was 
completely  coated  with  a  perfect  conductor.  Note  that  there  is  not  complete  spherical 
coverage  for  the  C-29  ISAR  data.  The  scattering  data  only  exists  for  limited  angular 
coverage,  and  this  is  depicted  in  Figure  2.4.  Also,  the  raw  ISAR  data  is  given  as 
samples  of  the  scattering  coefficients,  as  defined  in  Equation  (2.7). 

The  models  derived  from  the  raw  ISAR  data  did  not  utilize  all  of  the  data,  and 
the  parameters  of  the  data  used  to  develop  the  models  in  C29modeLm  are  shown  in 
Table  2.3.  It  was  assumed  for  all  the  C-29  models,  that  ten  scattering  centers  would 
be  used  to  model  the  scattering  data,  thus  F  =  10.  As  an  example,  lets  consider 
Model  Number  1.  For  Model  Number  1,  the  2-D  Damped  Exponential  Parame¬ 
ters,  determined  by  the  2-D  TLS-Prony  Technique  [3],  are  given  in  Equations  (2.11), 
(2.12),  and  (2.13).  In  these  equations,  the  angular  units  are  radians  in  the  expo¬ 
nential  arguments,  and  the  angles  range  from  — tt  to  tt.  Also,  in  the  MATLAB  code 
C29modeLm,  the  C-29  models  have  more  precision  than  five  digits  as  shown  in  the 
equations. 
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Figure  2.4:  Limits  of  angular  coverage  for  raw  C-29  ISAR  data. 


Table  2.2:  C-29  ISAR  Data  Specifications,  Unsealed. 


File  Name 

Freq  (GHz)  Range,  Step 

AZ  (Deg)  Range,  Step 

EL  (Deg) 

C29_3DJihN.02700 

26  to  36,  0.01 

-5  to  5,  0.04 

3 

C29_3DJihN.02500 

26  to  36,  0.01 

-5  to  5,  0.04 

5 

C29_3D_hhN.02300 

26  to  36,  0.01 

-5  to  5,  0.04 

7 

C29_3DJihW.02700 

26  to  36,  0.01 

19  to  29,  0.04 

3 

C29_3DJihW.02500 

26  to  36,  0.01 

19  to  29,  0.04 

5 

C29_3DJihW.02300 

26  to  36,  0.01 

19  to  29,  0.04 

7 

C29_3DJihW.02700 

26  to  36,  0.01 

-19  to  -29,  0.04 

3 

C29_3D_hhW.02500 

26  to  36,  0.01 

-19  to  -29,  0.04 

5 

C29_3DJihW.02300 

26  to  36,  0.01 

-19  to  -29,  0.04 

7 
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Table  2.3:  C29model.m  Data  Specifications,  Scaled 


Model  No. 

Freq  Range,  Steps  (GHz) 

Center  AZ,  Steps  (Deg) 

EL  (Deg) 

(M  X  N) 

1 

8.67  to  9.09,  0.0067 

-3.76,  0.04 

3 

64  X  64 

2 

8.67  to  9.09,  0.0067 

-0.08,  0.04 

3 

3 

8.67  to  9.09,  0.0067 

3.72,  0.04 

3 

4 

8.67  to  9.09,  0.0067 

-3.76,  0.04 

5 

64  X  64 

5 

8.67  to  9.09,  0.0067 

-0.08,  0.04 

5 

64  X  64 

6 

8.67  to  9.09,  0.0067 

3.72,  0.04 

5 

64  X  64 

7 

8.67  to  9.09,  0.0067 

-3.76,  0.04 

7 

64x64 

8 

8.67  to  9.09,  0.0067 

-0.08,  0.04 

7 

64  X  64 

9 

8.67  to  9.09,  0.0067 

3.72,  0.04 

7 

64  X  64 

10 

8.67  to  9.09,  0.0067 

20.24,  0.04 

3 

64  X  64 

11 

8.67  to  9.09,  0.0067 

23.92,  0.04 

3 

64  X  64 

12 

8.67  to  9.09,  0.0067 

27.72,  0.04 

3 

64  X  64 

13 

8.67  to  9.09,  0.0067 

20.24,  0.04 

5 

64  X  64 

14 

8.67  to  9.09,  0.0067 

23.92,  0.04 

5 

64  X  64 

15 

8.67  to  9.09,  0.0067 

27.72,  0.04 

5 

64  X  64 

16 

8.67  to  9.09,  0.0067 

20.24,  0.04 

7 

64  X  64 

17 

8.67  to  9.09,  0.0067 

23.92,  0.04 

7 

64  X  64 

18 

8.67  to  9.09,  0.0067 

27.72,  0.04 

7 

64x64 

19 

8.67  to  9.09,  0.0067 

-20.24,  0.04 

3 

20 

8.67  to  9.09,  0.0067 

-23.92,  0.04 

3 

64  X  64 

21 

8.67  to  9.09,  0.0067 

-27.72,  0.04 

3 

64  X  64 

22 

8.67  to  9.09,  0.0067 

-20.24,  0.04 

5 

64  X  64 

23 

8.67  to  9.09,  0.0067 

-23.92,  0.04 

5 

64  X  64 

24 

8.67  to  9.09,  0.0067 

-27.72,  0.04 

5 

64  X  64 

25 

8.67  to  9.09,  0.0067 

-20.24,  0.04 

7 

64  X  64 

26 

8.67  to  9.09,  0.0067 

-23.92,  0.04 

7 

64  X  64 

27 

8.67  to  9.09,  0.0067 

-27.72,  0.04 

7 

64  X  64 

16 


Px  = 


ai 

-0.1187+ O.lSllj 

0.1768e^(2-3067) 

a2 

0.0016  +  O.OOOOi 

0.0018e^(°-®^^^) 

as 

0.0001  -  O.OOOOj 

0.0001e^(-“-2738) 

a4 

0.1069 +  0.1440j 

0.1793e^'(°-^32i) 

as 

-0.0002 -0.1779i 

0.1779e^(-i-®'^^^) 

ae 

0.0008 +  0.0812j 

0.0812e^(^-5®^3) 

ar 

-0.0316  -  0.0104; 

0.0332e^'(-2-8244) 

as 

-0.0701  +  0.0726; 

0.10096^(2-3385) 

ag 

-0.0108  -  0.0978; 

0.0984e^("^-®^°®) 

aio 

-0.0694+0.0745; 

0.10186^(2.3207) 

Pxi 

1.0073-  0.1202; 

1.0145eJ(-°-^^^®) 

Px2 

0.3435  -  1.0489; 

1.10376^(-^-2344) 

PX3 

0.1962  -  1.1428; 

1.15956^(-i-'‘°°®^ 

Px4 

0.9399  -  0.3710; 

1.0104e^(-0-37®9) 

Px6 

0.9856  +  0.2002; 

1.00576^'(°-2004) 

PX0 

0.8361  -  0.5852; 

1.02056^'(-“-®i°^( 

PXT 

-0.1659-  1.0266; 

1.03996^(-i-'^310) 

Pxs 

1.0073-0.1202; 

1.01456^("°-^^®®) 

Px9 

0.9399  -  0.3710; 

1.01046^(-"-3^^®) 

Pxio 

-0.1206  +  0.9831; 

0.99056^(1-®®28) 

(2.12) 
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Pyi 

0.9246 -  0.3891i 

1.0031e(-°-3983) 

Py2 

1.0017-1- 0.0050j 

1.0017e(°‘^”^®) 

Pys 

1.0017-1-  0.0050i 

1.0017e(°“o^®) 

Py4 

1.0017 -f  0.0050j 

1.0017e(°““^®) 

Pys 

1.0017 -h0.0050i 

1.0017e(°“^^^ 

Pye 

0.9246-0.3891; 

1.0031e(-‘’-3^®3) 

Pyr 

1.0017-^0.0050; 

1.0017e<°’“4^) 

Pys 

1.0017 -f  0.0050; 

1.0017e(°'“^®) 

Pys 

0.9542  -  0.2781; 

0.9939e(-o-2336) 

Pyio 

0.9577-1-  0.3217; 

1.0103e(°-^2^^) 

(2.13) 


To  calculate  the  scattering  coefficient,  S  (/,  6)  =  d(m,  n),  at  EL  =  3°,  AZ  =  —4°, 
and  /  =  9GHz^  we  use  Equations  (2.3)  and  (2.11),  to  obtain 

m  =  50.0 
n  =  26.0 

and  thus 

S  {9GHz,  3°)  =  d(50, 26)  =  0.3577  +  0.3041 j  =  0.4695e^'‘^-^°^^.  (2.16) 

This  example  shows  how  to  calculate  the  scattering  coefficient  for  the  target  at  a 
given  frequency  and  aspect  angle.  The  angular  variable  9  used  here  directly  corre¬ 
sponds  to  AZ  angle  for  the  C-29.  To  get  a  better  appreciation  for  the  2-D  damped 
exponential  model,  consider  Model  Number  1  in  Table  2.3.  Using  the  raw  ISAR 
data  used  to  calculate  Model  Number  1  and  the  2-D  Inverse  Fast  Fourier  Transform 
(IFFT),  an  ISAR  image  was  formed  (see  Appendix  A.)  and  is  shown  in  Figure  2.5. 
Note  that  this  image  was  formed  from  actual  ISAR  data,  not  the  2-D  Damped  Ex¬ 
ponential  Model.  Also  note  that  an  outline  of  the  C-29  Aircraft  is  shown  on  the  plot 


(2.14) 

(2.16) 
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to  orient  the  reader  as  to  where  the  scattering  is  originating  from  on  the  C-29  scale 
model.  The  dimensions  of  the  image  plots  have  been  scaled  to  full-scale  C-29  size. 

Next,  the  2-D  Damped  Exponential  Model,  Model  Number  1,  was  used  to  generate 
ISAR  data  over  the  same  angular  and  frequency  ranges  as  C-29  Model  Number  1,  and 
a  model  based  image  was  formed.  This  image  is  shown  in  Figure  2.6.  As  before,  an 
overlay  of  the  C-29  is  shown.  Also,  the  scattering  center  locations  of  the  ten  scattering 
centers  defined  by  the  angles  of  Px  and  py,  as  given  in  Equations  (2.12)  and  (2.13),  are 
overlayed  on  the  two  images  and  denoted  by  ’o’s.  Comparing  Figures  2.5  and  2.6,  it 
can  be  seen  the  2-D  Damped  Exponential  Model  does  not  reproduce  an  exact  image, 
but  does  reproduce  an  image  which  is  visually  close  to  the  actual  ISAR  image.  Also, 
the  scattering  centers  are  placed  on  the  ISAR  image  at  locations  where  most  of  the 
scattering  energy  is  being  reflected  from. 

Using  the  2-D  Damped  Exponential  (2-D  Prony’s)  Model,  many  different  radar 
targets  can  be  simulated  by  placing  scattering  centers  in  various  locations.  The  MAT- 
LAB  code  point_inodel.m  creates  a  generic  aircraft  by  placing  scattering  centers  at 
probable  scattering  locations.  This  ’’pointjnodel”  can  be  used  by  profile. m  to  gen¬ 
erate  an  entire  radar-target  engagement.  Now  that  the  scattering  characteristic  of  the 
target  have  been  described,  the  scattered  field  can  easily  be  calculated,  as  discussed 
in  the  next  section. 

2.4  Scattered  Field  Calculations 

To  calculate  the  scattered  field  at  time  t,  we  can  use  the  variables  in  Table  2.1 
along  with  Equations  (2.1)  and  (2.4).  First,  define  the  distance  from  the  target  phase 
center  to  the  radar  receive  antenna  as 

r  =  \Rt-  Ro\  (2.17) 

where  |  •  |  stands  for  spatial  distance  in  this  case.  Using  the  definitions  of  Ptr  and 
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Actual  SAR  Image,  Radar  at  Top  Ctr,  AZ=-3.762,  EL=3.007 


-10  -5  0  5  10 

Cross-Range  Dimension,  meters 


Figure  2.5;  Actual  ISAR  Image  for  One-Third  C-29  Scale  Model.  An  overlay  of  the 
C-29  Aircraft  is  shown.  Dimensions  have  been  scaled  to  correspond  to  a  full-size  C-29 
Aircraft. 
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Down-Range  Dimension,  meters 


Model  Gen  SAR  Image,  Radar  at  Top  Ctr,  AZ=-3.762,  EL=3.007 
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Figure  2.6:  Model  Generated  ISAR  Image  for  One-Third  C-29  Scale  Model.  An 
overlay  of  the  C-29  Aircraft  is  shown.  Dimensions  have  been  scaled  to  correspond  to 
a  full  size  C-29  Aircraft. 
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Gtr,  the  complex  frequency  domain  value  of  the  scattered  field  as  seen  by  the  receiver 
is 


®scat  —  n) 


f  2{Ptr){Gtr)r]o 
\  dTrr^ 


,-jko2r 


(2.18) 


where 


K  =  ^  (2.19) 

c 

rjo  =  376.7  n 
c  =  2.9979  X  10®  m/s. 


This  equation  was  derived  assuming  a  mono-static  radar,  and  that  the  velocity  of  the 
aircraft  was  much  smaller  than  the  speed  of  light  (which  is  a  pretty  good  assumption). 
From  this  complex  quantity,  an  actual  time  domain  expression  for  the  scattered  field 
can  be  found.  Let 

Escat  =  |Escat|e^’^®““*,  (2.20) 

where  |  •  |  denotes  complex  magnitude  and  L  denotes  complex  angle.  Next,  the 
dependence  is  added  back  into  the  frequency  domain  expression  to  yield 


EL,  =  E,eate^“*. 

Next,  define 

Eo  —  lEgcatl 

4^0  —  7Eg,,j,, 


(2.21) 


(2.22) 

(2.23) 


where  Eo  and  (po  are  the  magnitude  and  phase  of  the  complex  number  Finally, 

a  time  domain  expression  for  the  scattered  field  is  given  by 


Eg{t)  -  Eo  cos  {ut  -h  ZEscat)  =  Eo  cos  {(po)-  (2.24) 
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This  time  domain  expression  is  valid  for  time  t,  frequency  f  {co  =  27r/),  and  the 
physical  situation  described  above.  After  the  target  and/or  the  radar  move,  Egcat 
will  change  its  value  and  thus  change  the  values  and  Es{t). 

Some  important  assumptions  have  been  made  in  order  to  calculate  the  scattered 
field  in  Equation  (2.24).  First,  the  2-D  Damped  Exponential  Model  is  physically,  by 
definition,  a  two-dimensional  model.  Thus,  the  aircraft  is  modeled  in  space  by  a  finite 
plate,  or  plane  in  3-D  space.  The  direction  of  this  plane  is  defined  by  the  directions 
Pt  and  a  line  parallel  to  the  s2:-plane,  as  demonstrated  in  Figure  2.1.  For  the  simple 
case  of  Pt  lying  in  the  (— a:)-direction,  the  aircraft  model  is  parallel  to  the  a; plane. 
Second,  all  of  the  scattering  centers  lie  in  this  plane  within  the  physical  boundaries 
of  the  unambiguous  range  of  the  target.  The  unambiguous  range  is  defined  by  the 
frequency  sampling  rate  of  the  raw  ISAR  data,  df  and  da.  The  unambiguous  range 


in  the  down-range  dimension  is 


Ruf  = 


2(d/) 


(2.25) 


while  the  unambiguous  range  in  the  cross-range  dimension  of  the  image  is 


(2.26) 


where  dfa  =  {da)radf  {darad  is  the  value  of  da  in  radians).  The  down-range  and  cross- 
range  unambiguous  ranges  are  determined  by  Nyquist’s  sampling  criteria  on  the  raw 
ISAR  data  in  the  2-D  frequency  domain.  Also  note  that  in  the  ISAR  images  shown  in 
Figure  2.5,  the  unambiguous  ranges  are  the  limits  on  the  horizontal  (cross-range)  and 
vertical  (down-range)  axes.  If  the  actual  target  size  exceeds  the  unambiguous  ranges, 
aliasing  will  occur  in  the  ISAR  image  and  2-D  Damped-Exponential  parameters  can 
not  be  estimated. 


Now  that  the  scattering  characteristics  of  the  target  are  known,  aircraft  motions 
and  vibrations  can  be  added  to  the  nominal  flight  path  to  better  simulate  a  dynamic 
aircraft  in  flight. 


2.5  Aircraft  Vibrations  and  Motion  Models 


It  is  weU  know  that  any  aircraft  vibrates  slightly  during  flight.  Actually,  there  are 
many  modes  of  vibrations.  These  vibrations  are  significant  factors  with  respect  to 
the  Doppler  shift  and  angular  error  as  seen  by  a  radar  tracking  a  target.  There  are 
two  types  of  aircraft  motion  simulated  here.  The  first  is  called  rigid-body  motion, 
and  the  second  is  called  individual  motion.  Rigid-body  motion  refers  to  the  entire 
aircraft  vibrating,  or  moving,  away  from  the  intended  path  of  flight  as  a  rigid-body. 
In  this  case  there  is  no  flexing  of  the  aircraft  at  all.  The  second  case  involves  flexing 
of  the  aircraft.  In  this  case  the  aircraft  flexes  and  the  various  parts  of  the  aircraft 
move  with  respect  to  one  another.  Next,  each  type  of  motion  is  described. 

2.5.1  Rigid  Body  Motion 

In  rigid  body  motion,  there  is  no  flexing  of  the  aircraft  structure  at  all.  There 
are  two  types  of  rigid  body  motion,  translational  and  rotational.  For  translational 
motion,  the  entire  aircraft  translates  from  its  nominal  flight  path  in  all  three  Cartesian 
directions  at  each  instant  of  time.  For  rotational  motion,  the  entire  aircraft  rotates 
from  its  direction  in  all  three  angular  directions  (roll,  pitch,  and  yaw).  The  MATLAB 
code  which  calculates  the  translational  and  rotational  rigid-body  motion  is  called 
rigid_motion.m.  First,  lets  consider  translational  motion. 

Translational  Rigid  Body  Motion 

The  aircraft  position  with  no  motion  is  given  by  Rt  for  a  given  time  t.  Transla¬ 
tional  motion  changes  the  (a;,  y,  z)  components  of  Rt  by  slight  amounts  to  simulate  the 
motion.  Consider  three  independent  identically  distributed  white  Gaussian  random 
processes,  wtny,  and  wtn^,  each  with  zero  mean  and  standard  deviation  cr^. 

Assuming  there  are  ti  time  steps  in  the  simulation,  each  process  has  ti  samples.  First, 
each  process  is  filtered  by  a  16-th  order  FIR  low-pass  filter  with  cut-off  frequency  ft^ 
with  discrete-time  impulse  response  htr[n].  The  resulting  three  processes  are  ctrix, 


24 


ctfiy^  Biiid.  ctTi^^  defined  by 


CtUx  —  WtUx  0  hfr 

(2.27) 

CtUy  =  WtUy  0  htT 

(2.28) 

ctrix  =  wtn^  0  htr 

(2.29) 

where  0  denotes  discrete-time  convolution  [13].  Next  a  3  x  3  correlation  matrix  is 
defined  as 


pxx  pxy  pxz 
Pyx  Pyy  pyz 


pxx  Pzy  Pzz  J 

At  each  instant  of  time,  three  new  random  variables  are  defined  as 


(2.30) 


Ctx  [^] 

ctnx[n] 

cty[n] 

=  C', 

ctny[n\ 

ct,\n] 

ctn^[n] 

(2.31) 


Thus,  the  three  new  random  processes,  c4,  cty,  and  ct,  are  colored  (by  filtering)  and 
correlated  (by  Equation  (2.31))  with  one  another.  At  time  t,  the  new  aircraft  position 


Bit  —  Rt-{- 


ct^  [n] 

cty  [ra] 

ctz[n] 


(2.32) 


Thus,  the  aircraft  is  translationaUy  shifted  from  the  intended  flight  path  by  bandlim- 
ited  correlated  white  noise.  Note  that  the  magnitude  of  the  shifts  tend  to  be  smaller 
than  tr,*  since  the  low-pass  filtering  filters  out  the  high  frequency,  or  large  magnitude 


components.  The  cut-off  frequency,  ft^,  is  nominally  around  2  to  5  Hz.  The  standard 


deviation  of  the  translational  movement,  is  around  10  cm  per  dimension,  while 
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the  correlation  matrix  is  typically 

0.9  0.1  0.1 

C'i=  0.1  0.9  0.1  .  (2.33) 

0.1  0.1  0.9 

All  of  these  values  were  chosen  by  consulting  with  an  aeronautical  engineer  [14]. 
Rotation2il  Rigid  Body  Motion 

Next,  the  rotational  rigid  body  motion  is  defined.  This  consists  to  perturbations 
the  physical  direction  the  aircraft  is  pointing.  Recall  the  unit  vector  Dt  is  the  direction 
the  aircraft  is  traveling  in  while  the  unit  vector  Pt  is  the  direction  the  aircraft  is 
physically  pointing  in.  See  Figure  2.2  for  a  visual  definition  of  Pt.  Next,  define  the 
three  orthogonal  angular  components,  roll  IZ,  pitch  P,  and  yaw  ^  as  the  standard 
roll,  pitch,  and  yaw  angles  associated  with  an  aircraft  in  flight.  Defining  the  angular 
vector  A, 

n 

V  .  (2.34) 

y 

If  ^  =  [0  0  0]^,  where  T  denotes  transpose,  then  roll,  pitch,  and  yaw  are  all  zero. 
In  this  case,  the  aircraft  centerline  is  pointing  in  the  Dt  direction  and  Dt  =  Pt.  When 
roll,  pitch,  and  yaw  are  not  zero,  then  .4  [0  0  0]^,  and  Dt  ^  Pt.  In  this  case,  the 

unit  vector  Pt  is  pointing  in  the  direction  defined  by  the  roll,  pitch,  and  yaw  values. 
The  aircraft  is  always  traveling  in  the  Dt  direction,  but  is  always  pointing  in  the  Pi 
direction. 

Thus,  Pt  is  a  function  of  Dt  and  A.  When  Dt  ^  Pt,  the  aircraft  is  no  longer 
oriented  in  the  Dt  direction,  and  new  azimuth  (AZ)  and  elevation  (EL)  angles  need 
to  be  found  (see  Figure  2.2).  The  MATLAB  code  scattered-field. m  calculates  the 
new  AZ  and  EL  angles  due  to  any  non-zero  values  in  A.  The  equations  for  the  new 
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AZ  and  EL  angles  are  not  given  here,  but  are  given  in  scattered-field. m,  shown  in 
Appendix  B. 

As  with  translation  motion,  rotational  rigid-body  motion  consists  of  adding  cor¬ 
related  bandlimited  noise  to  the  roll,  pitch,  and  yaw,  or  to  Pt.  Since  the  movement 
values  are  low-pass  filtered,  the  actual  angular  perturbations  tend  to  be  smaller  than 
Urri  the  standard  deviation  of  the  original  white  Gaussian  processes  for  rotational  rigid 
body  motion.  Defining  the  movement  vector  in  the  same  fashion  as  for  translational 
motion. 


cfj,[n] 
cf\n\  =  I  cry\n] 
cr^[n\ 

and  the  new  aircraft  pointing  direction  becomes  (at  time  step  n) 


(2.35) 


Pt  =  Pt  +  cf[n]. 


(2.36) 


The  standard  deviation  of  the  rotational  movement,  cr^r,  is  typically  around  2° 
and  bandlimited  to  2  to  5  Hz  in  each  angular  dimension.  The  correlation  matrix  for 
rotational  motion  is  identical  to  the  one  for  translational  motion  [14]. 

Translational  and  rotational  rigid  body  motion  are  critical  to  the  dynamic  simula¬ 
tion  of  a  radar-target  encounter  since  this  motion  has  a  significant  effect  on  scattered 
field  values,  Doppler  shift,  and  angular  error.  The  other  type  of  motion  considered 
is  called  individual  scattering  center  motion.  Unlike  rigid-body  motion,  in  individual 
scattering  center  motion  the  scattering  centers  on  the  target  move  independent  to 
one  another.  This  type  of  motion  is  described  next. 


2.5.2  Individual  Scattering  Center  Motion 

In  order  to  simulate  flexing  of  the  aircraft  structure,  which  does  occur,  individual 
movement  of  the  scattering  centers  is  implemented  in  the  simulations.  The  scatter¬ 
ing  centers  are  allowed  to  move  in  the  physical  2-D  plane  in  which  the  scattering 
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centers  reside.  Thus,  the  down-range  and  cross-range  locations  of  the  scattering 
centers  are  varied.  Assuming  that  there  are  F  scattering  centers,  2V  discrete  time, 
independent  identically  distributed  Gaussian  random  processes  are  created,  each  of 
length  ti.  The  standard  deviation  of  all  of  these  processes  is  aim-  The  processes 
are  called  and  {uIMxri}J^-^.  Next,  each  of  the  random  processes  are 

filtered  through  a  16-th  order  discrete-time  low-pass  FIR  filter  to  produce  the  colored 
sequences  and  The  cut-off  frequency  for  the  filter  is  fim^. 

Next,  it  is  desired  to  change  the  positions  of  the  individual  scattering  centers  by 
adding  the  values  of  IMdvi  and  IMxVi  to  the  positions  of  the  individual  scattering 
centers.  Recall  from  Section  2.2  that  the  angles  of  the  x  and  y-poles  determine  the 


locations  of  the  scattering  centers  in  the  ISAR  image.  Thus,  the  angles  of  the  x  and 
y-poles  will  be  individually  changed  in  order  to  move  the  scattering  centers  around  on 
the  image.  Also,  for  a  given  instant  in  time  (n)  and  one  scattering  center,  the  down- 
range  and  cross-range  movement  distances  are  given  by  IMdri[n]  and  IMxri[n\. 
However,  before  these  perturbations  are  applied,  the  down-range  and  cross-range 
components  are  correlated  using  a  2  x  2  correlation  matrix,  yielding  cIMdvi  and 
cIMxTi.  Thus, 


cIMdri 

cIMxVi 


Pdd  Pdx 
Pxd  Pxx 


IMdvi 

IMxvi 


(2.37) 


The  movement  of  the  individual  scattering  centers  will  change  the  scattered  field  from 
the  target.  Recall  that  the  scattered  field  is  calculated  using  Equation  (2.18).  The 
value  of  d(m,  n)  will  change  if  the  scattering  centers  move  since  movement  of  the 
scattering  centers  is  simulated  by  changing  the  angles  of  the  x  and  ?/-poles. 

As  an  example,  consider  50  overlayed  realizations  of  individual  scattering  center 
movement  as  seen  on  the  ISAR  images  shown  in  Figures  2.7  and  2.8.  This  corresponds 
to  how  the  individual  scattering  centers  move  as  the  simulation  progresses  through 
50  time  increments.  In  Figure  2.7,  the  standard  deviation  of  movement,  is  set  to 
2  cm.  The  scattering  center  locations  are  denoted  by  ’-f’s  while  the  original  locations 


28 


are  denoted  by  ’o’s.  Figure  2.8  shows  the  case  when  (7^^,  ==  20cm.  For  both  cases, 
firric  =  20Hz  and  the  correlation  matrix  is  given  by 


Pdd  Pdx 
Pxd  Pxx 


0.8  0.2 

0.2  0.8 


(2.38) 


In  each  figure  the  model  based  ISAR  image  is  shown,  generated  from  the  “un-moved” 
2-D  Damped  Exponential  parameters.  An  aircraft  overlay  is  also  shown  with  each 
plot. 

Please  note  that  20cm  is  somewhat  large  for  the  magnitude  of  the  individual  mo¬ 
tion,  but  it  is  shown  here  to  demonstrate  how  the  scattering  centers  change  position. 
More  realistic  movement  numbers  are  1  to  5cm,  values  within  this  range  are  used  for 
the  simulations  in  Section  4. 

The  individual  scattering  center  motion  developed  here  is  implemented  in  the 
MATLAB  function  indiv_mot.m.  This  type  of  motion  is  critical  to  achieve  Doppler 
signatures  and  angulax  errors  which  are  realistic,  as  shown  in  Section  4.  It  is  important 
to  note  that  the  individual  scattering  center  motion  values  do  not  depend  upon  where 
on  the  target  the  scattering  center  is  located.  Thus,  a  scattering  center  at  the  phase 
center  of  the  target  will,  statistically,  be  moved  just  like  a  scattering  center  near  the 
edge  of  the  unambiguous  range.  As  an  upgrade  to  this  simulation  package,  some 
weighting  factor  could  be  imposed  upon  the  magnitude  and  frequency  of  individual 
scattering  center  motion  as  a  function  of  each  scattering  center’s  location.  This  is 
discussed  more  in  Section  6. 


2.6  Summary  of  Dynamics  and  Target  Scattering 

In  summary,  the  typical  radar-target  engagement  consists  of  a  stationary  radar 
(which  could  be  a  missile  seeker)  positioned  at  the  origin  in  Figure  2.1  and  a  target 
(C-29  aircraft)  flying  in  3-D  space.  As  the  aircraft  is  flying  straight  and  level,  it  is 
vibrating  both  as  a  rigid-body  and  as  a  flexing  aircraft.  At  discrete  time  intervals,  dt, 
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Model  Gen  SAR  Image,  Radar  at  Top  Ctr,  AZ=-3.762,  EL=3.007 
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Figure  2.7:  Individual  Scattering  Center  Motion  Plots.  Plot  shows  2cm  motion  with 
50  overlayed  realizations.  An  aircraft  overlay  is  shown  with  the  plot. 
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Model  Gen  SAR  Image,  Radar  at  Top  Ctr,  AZ=-3.762,  EL=3.007 
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Figure  2.8:  Individual  Scattering  Center  Motion  Plots.  Plot  shows  20cni  motion  with 
50  overlayed  realizations.  An  aircraft  overlay  is  shown  with  the  plot. 
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the  simulation  calculates  the  scattered  electric  field,  the  Doppler  frequency  shift,  and 
the  angular  error  (glint)  all  as  seen  at  the  radar.  In  the  next,  all  of  the  simulation 
outputs,  scattered  field,  Doppler  shift,  and  angular  error  are  described  in  detail. 


3.  Simulation  Outputs 

In  this  section,  the  critical  output  quantities  of  the  simulations  are  derived  and 
defined.  They  consist  of  the  scattered  electric  field,  the  Doppler  shift,  and  the  angular 
error  (glint)  as  seen  at  the  radar  receiver.  Each  of  these  quantities  is  calculated  at 
each  time  step  in  the  simulation. 

3.1  Scattered  Field 

The  scattered  field  is  a  critical  output  quantity  as  this  is  the  value  of  the  electric 
field  as  seen  at  the  radar  receiver.  In  a  hardware  in  the  loop  simulation,  this  number 
gives  the  value  of  the  electric  field  which  should  be  induced  at  the  receiver  at  time 
t.  The  polarization  of  the  electric  field  will  be  defined  by  the  polarization  of  the 
scattering  coefficient  in  Equation  (2.2).  For  the  C-29  hh-polarization  data,  the  scat¬ 
tered  field  (and  the  incident  field)  is  horizontally  polarized.  The  actual  value  of  the 
scattered  field  was  developed  in  Section  2.4,  and  is  given  by  Es{t)  in  Equation  (2.24). 
This  value  is  calculated  in  the  MATLAB  code  scatteredJield.m  and  is  passed  to 
profile.m. 

If  the  magnitude  and  phase  of  the  scattered  field  are  desired  as  a  function  of  time, 
then  they  are  also  given  in  Section  2.4  as  Eo  and  in  Equations  (2.23)  and  (2.23),  re¬ 
spectively.  These  quantities  are  calculated  in  the  MATLAB  code  scattered_field.m 
but  are  not  passed  to  profile.m.  They  are  not  part  of  the  standard  output  of  the 
simulation  package,  but  if  desired  can  be  extracted  from  scattered_field.m. 
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Another  useful  quantity  related  to  the  scattered  electric  field  is  the  power  received 
at  the  radar,  denoted  Pr.  Assuming  the  radar  receive  antenna  gain  is  the  same  as 
the  transmitter  antenna  gain,  which  is  appropriate  for  a  mono-static  radar,  the  power 
received  immediately  after  the  radar  receive  antenna  is  given  by 

=  (3.1) 

where  |  •  |  denotes  absolute  value,  rjo  and  c  are  defined  in  Equation  (2.20),  Gtr  is 
defined  in  Table  2.1,  and  /  is  the  frequency.  This  value  is  also  calculated  in  the 
MATLAB  code  scatteredJield.ni. 

3.2  Doppler  Shift 

The  Doppler  shift  is  an  important  quantity  for  radar  receivers  since  this  will  define 
the  “Doppler  Bandwidth”  over  which  a  radar  must  operate.  The  Doppler  shift  refers 
to  the  instantaneous  frequency  of  the  scattered  electric  field  at  the  radar  receiver. 
The  standard  definition  of  Doppler  shift  is  the  difference  between  the  frequency  of 
the  electric  field  with  radar-target  motion  and  that  without  radar-target  motion.  For 
example,  if  a  CW  radar  is  illuminating  a  target  with  an  /  =  9  GHz  electric  field,  and 
the  target  is  moving  at  100  meters/sec  directly  toward  the  radar,  then  the  Doppler 
shift  would  be  given  by  [4] 

/i  =  +^^^^  =  +6  kHz.  (3.2) 

c 

The  frequency  of  the  scattered  field  received  by  the  radar  is  given  by 

fr  =  f  +  fd-  (3.3) 

For  this  numerical  example,  the  received  frequency  would  be  9.000006  GHz.  This 

scenario  is  a  very  simple  case.  When  both  the  target  and  radar  are  moving,  along 
with  different  parts  of  the  target  moving  at  different  rates  (individual  motion),  the 
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calculation  of  the  Doppler  shift  is  a  bit  more  involved.  We  will  derive  a  general 
expression  for  the  Doppler  shift  as  follows. 

Starting  with  Equation  (2.18),  we  can  rewrite  the  frequency  domain  expression 
for  the  scattered  electric  field  for  a  given  frequency  /  (normalized  to  m)  and  aspect 
angle  6  (normalized  to  n)  at  time  t  as 


(3.4) 

(3.5) 


r  =  number  of  scattering  centers  (3.6) 

Ei  =  amplitude  of  i-th  scattering  center  (3.7) 

(/•i  =  phase  of  the  i-th  scattering  center  (3-8) 

rii  =  unit  vector  pointing  from  z-th  scattering  center  to  radar  (3.9) 

Rot  =  Rt  —  vector  from  phase  center  of  target  to  radar  (3.10) 

a;  =  radian  frequency(c<;  =  27r/)  (3-11) 

•  =  standard  vector  dot  product.  (3.12) 


It  is  advantageous  to  rewrite  the  frequency  domain  expression  for  the  scattered  field 
in  this  form.  It  must  be  emphasized  the  the  expression  for  Eg^at  in  Equation  (3.5)  and 
the  expression  for  in  Equation  (3.4)  yield  the  same  value,  they  are  just  different 
methods  of  expressing  the  same  quantity.  Also,  the  values  for  in  both  of  those 
equations  only  differ  from  the  value  for  Escat  in  Equation  (2.18)  by  an  term.  The 
factor  of  2  multiplying  Rot  accounts  for  the  two-way  travel  of  the  electromagnetic 
fields  (radar  to  target  and  target  to  radar).  Note  that  now  the  scattered  field  has 
been  expressed  as  the  sum  of  the  scattered  fields  from  each  of  the  T  scattering  centers. 
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This  representation  is  useful  for  developing  an  expression  for  Doppler  shift  and  also 
for  developing  an  expression  for  angular  error,  as  will  be  shown  in  Section  3.3.  Next, 
bringing  the  term  outside  of  the  summation  we  have 

ELat  =  (3.13) 

i~l 

Consider  the  phase  angle  of  and  define 

^  -  cut  ^  (f>  -  ojt  =  ZEscat  (3.14) 

where  (j)  was  defined  in  Equation  (2.23).  Thus,  $  is  the  phase  of  the  complex  number 
Escat-  The  Doppler  shift,  fd,  is  now  defined  as 

(3-15) 

When  there  is  no  target-radar  motion,  the  time  derivative  of  ^  is  zero  and  so  is  fd- 
However,  when  there  is  radar-target  motion,  ^  is  not  equal  to  zero  and  we  have  a 
Doppler  shift. 

To  calculate  let 

ELi  =  e’"‘  (3.16) 

where  Eo  (defined  in  Equation  (2.23))  and  $  are  real  numbers  and  are  the  complex 
magnitude  and  phase,  respectively,  of  with  the  e-'"*  dependence  removed.  Next, 
let 

Eoe^^  =  G  +  jF  (3.17) 

where  G  and  F  are  both  real  numbers.  Further,  F  and  G  can  be  written  as 

G  =  REAbjEoe^*}  =  Ef=iF^iCos(0i  -  2A:<,nJ  •  ifot) 

F  =  IMAG|Eoe-^*|  =  sin.  {(f)i  -  2ko^i  ■  Rot) 
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where  the  REAL  operator  takes  the  real  part  of  a  complex  number  and  the  IMAG 
operator  takes  the  imaginary  part.  Then 


#  =  arctan 


(3.19) 


Taking  the  time  derivative  of  $ 


dt 


d  r 

—  <  arctan 

dt  I 


GF'  -  FG' 

G2 


(3.20) 


where  F'  is  the  time  derivative  of  F  and  G'  is  the  time  derivative  of  G.  Rearranging 
terms  and  realizing  that 


5=  1  (GF'-FG'). 


(3.21) 


In  the  most  general  case,  assuming  that  Ui,  and  Rot  are  all  functions  of  time,  the 
time  derivative  of  $  reduces  to 


1^1  d(j)i 


-2ko 


dui 


dt  Eo^  ^  *  1  °  \  dt 


Rot  + 


dR, 


oi 


dt 


m 


(3.22) 


where 


i 

ai  =  EiY,  cos  (j>i  -  2ko  (n^  •  Rot)  -  <t>i  +  2ko  {ni  •  Ro?j)  ■  (3.23) 


1=1 


Finally,  using  Equations  (3.15)  and  (3.22),  an  expression  for  the  Doppler  shift  is 


~  -p 

i=l 


d(f)i 

dt 


2ko 


drii 

dt 


Rot  + 


dR. 


ot 


dt 


Ui 


(3.24) 


The  Doppler  shift  is  calculated  by  Equation  (3.24),  where  at  is  defined  in  Equa¬ 
tion  (3.23),  in  the  MATLAB  code  scattered-field. m.  The  three  required  time  deriva¬ 
tives,  and  are  calculated  by  taking  the  difference  between  their  present 

values  and  their  values  at  the  previous  iteration,  and  dividing  by  the  time  step,  dt. 
This  is  the  most  accurate  method  available  for  calculating  the  Doppler  shift.  The 
example  in  Section  4  shows  how  the  Doppler  shift  varies  over  a  flight  profile.  The 
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range  over  which  the  Doppler  shift  varies  as  a  function  of  time  determines  the  required 
“Doppler  bandwidth”  of  a  radar  (or  missile)  receiver.  Also,  using  Equation  (3.3),  the 
frequency  of  the  scattered  field  received  by  the  radar,  f^,  can  be  determined. 


3,3  Angular  Error 


The  angular  error,  or  glint,  is  another  important  quantity  for  a  radar  or  missile 
receiver.  The  angular  error,  or  glint,  is  defined  as  the  angular  difference  between  the 
vector  from  the  target  phase  center  to  the  radar  receiver  and  the  direction  of  arrival 
of  the  electric  field  at  the  radar  receiver.  The  direction  of  arrival  (DOA)  is  defined  as 
the  direction  of  the  vector  normal  to  the  incoming  wavefront  of  the  scattered  electric 
field  [15,  16,  17].  The  DOA  can  also  be  defined  as  the  direction  of  the  Poynting  Vector 
of  the  Transverse  Electromagnetic  Wave  (TEM)  incident  at  the  radar  receiver  [17]. 
First,  an  expression  for  the  DOA  vector  will  be  found.  The  DOA  vector  is  a  vector 
pointing  in  the  DOA,  and  is  called  d. 

The  DOA  vector,  d,  is  related  to  the  gradient  of  the  phase  of  the  scattered  elec¬ 
tric  field  incident  at  the  radar  receiver.  Equation  (3.19)  gives  an  expression  for  the 
required  phase,  $.  Rewriting  this  expression  using  the  definitions  of  F  and  G,  we 
obtain 


$  =  arctan 


T,f:=iEi  sin  {(t>i  -  2koni  •  Rot) 


(3.25) 


COS  ((^j  2Ajo7lj  *  Rot)  ^ 

Using  Equations  (3.5)  and  (3.25),  the  DOA  vector  dis  related  to  the  phase  $  by  the 
gradient  function,  namely. 


r  ^0$ 


(3.26) 


where  V  is  the  gradient  operator  and  x,  y,  and  z  are  the  unit  vectors  in  the  x,  y,  and 
2  spatial  Cartesian  directions,  respectively.  The  partial  derivatives  of  $  with  respect 
to  the  X,  y,  and  z-directions  are  easily  computed.  For  a  detailed  computation  of  the 
partial  derivatives,  see  [16].  Taking  all  three  partial  derivatives  and  combining  the 
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results,  a  final  expression  for  the  DOA  vector  is 

=  (3.27) 

where  ctj  is  defined  in  Equation  (3.23)  and  nj  is  defined  in  Equation  (3.12).  The  factor 
of  2  in  front  of  kg  is  due  to  the  two-way  propagation  of  the  electric  field  from  radar 
to  target. 

Next,  recall  the  radar-target  direction  vector,  Rgt,  defined  in  Equation  (3.10),  as 
the  vector  from  the  radar  receiver  to  the  phase  center  of  the  target.  The  angular 
error,  Qg,  is  defined  as  the  angular  difference  between  the  vectors  d  and  Rot-  This 
angular  error  is  calculated  using  the  dot  product  as 

Q,e  =  arccos  (  C  ^  ]  (3.28) 

vMPotiy 

where  •  denotes  the  vector  dot  product  and  |  •  |  denotes  the  magnitude  of  the  vector. 
The  sign  of  Og  is  determined  using  the  cross  product  between  d  and  Rgt- 

The  angular  error  is  an  important  quantity  since  radar  receivers  and  missile  seekers 
point  in  the  direction  defined  by  the  gradient  of  the  incident  field,  the  direction  of  d, 
not  in  the  direction  of  the  target,  Rgt.  Note  that  Qg  is  typically  a  very  small  number. 
In  the  simulations,  fig  is  calculated  in  the  MATLAB  code  scatteredJield.m.  In  the 
simulations,  Qg  is  compared  with  the  apparent  angular  wingspan  of  the  C-29  aircraft 
as  seen  at  the  radar  receive  antenna. 


4.  Simulation  Example 

In  this  section,  a  typical  simulation  was  performed  and  the  output  quantities  are 
shown.  First,  the  profile  is  described.  Figure  4.1  shows  a  MATLAB  plot  of  the 
simulated  flight  path  of  the  C-29  aircraft.  The  radar  (or  missile  seeker)  is  located  at 
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Table  4.1:  Parameters  for  sample  simulation. 


Variable 

Description 

Value 

/ 

radar  frequency 

8.6674  GHz 

Ptr 

transmitter  power 

IHIIIIIQQ^HIII 

Gtr 

transmitter  antenna  gain 

Vt 

target  velocity 

150  m/s 

Dt 

target  direction  vector  (x,y,z) 

(1.0.0) 

Rti 

init.  target  position  (x,y,z) 

(2800,-100,-100)  m 

Ro 

radar  position  (x,y,z) 

(0,0,0)  m 

dt 

time  step  for  simulations 

0.01  s 

ti 

number  of  iterations 

601 

t 

time 

0  to  6  s 

n 

normalized  time 

0  to  601 

the  origin  of  the  coordinate  system.  Table  4.1  list  the  simulation  parameters.  Note 
that  the  values  shown  correspond  to  a  full-scale  C-29  aircraft.  The  aircraft  had  both 
rigid-body  motion  and  individual  scattering  center  motion.  For  the  translational 
rigid-body  motion,  the  standard  deviation  of  each  of  the  translational  components 
was  set  to  10  cm,  the  cut-off  frequency  to  2  Hz,  and  the  correlation  matrix  is  as  given 
in  Equation  (2.33).  For  the  rotational  rigid-body  motion,  the  standard  deviation  of 
each  of  the  angular  components  was  set  to  2°,  the  cut-off  frequency  to  2  Hz,  and  the 
correlation  matrix  is  as  given  in  Equation  (2.33).  The  individual  scattering  center 
motion  had  a  magnitude  of  2  cm  in  each  direction,  with  a  cut-off  frequency  of  10  Hz, 
and  the  down-range  and  cross-range  motion  components  are  correlated  by  the  matrix 
in  Equation  (2.38).  These  values  were  chosen  in  consultation  with  an  aeronautical 
engineer  [14]. 

The  results  of  the  simulation  are  shown  in  Figures  4.2  through  4.6.  Examining 
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Figure  4.2,  we  can  see  the  effects  of  the  C-29  vibrations  in  the  azimuth  and  elevation 
plots.  The  magnitude  of  the  angular  vibrations  is  smaller  than  the  Urr  =  2°  due 
to  the  low  pass  filtering  described  in  Section  2.5.1.  Since  the  radar-target  range  is 
much  greater  than  the  vibrations,  the  effects  of  the  vibrations  can  not  be  seen  in  the 
radar-target  range  plot.  Also  note  C-29  Model  Number  1,  as  defined  in  Table  2.3, 
was  the  best  model  to  use  for  this  simulation,  and  was  used  exclusively  throughout 
this  simulation  run.  Recall  that  the  center  AZ  and  EL  values  for  Model  Number  1 
are  —3.76°  and  3°,  respectively.  These  values  are  closer  to  the  required  AZ  and  EL 
values  shown  in  Figure  4.2  than  the  values  for  any  other  available  C-29  models  listed 
in  Table  2.3. 

Figure  4.3  shows  both  the  scattered  electric  field  at  the  radar  receiver,  Es{t)  as 
defined  in  Equation  (2.24),  plotted  in  dB,  and  the  angular  error  at  the  radar  receiver, 
fie  as  defined  in  Equation  (3.28),  plotted  in  degrees,  both  versus  time.  The  rapid 
variations  in  the  scattered  field  can  mainly  be  attributed  to  the  cos  [wt)  dependence 
of  Es{t).  The  other  variations  in  Es{t)  can  be  attributed  to  the  dynamics  of  the 
scenario,  along  with  the  rapidly  changing  scattering  Characteristics  of  the  C-29  as 
a  function  of  AZ  and  EL  angles.  The  angular  error,  fig?  shown  in  Figure  4.3,  is 
referenced  against  the  angular  extent  (or  wingspan)  of  the  C-29  aircraft  as  seen  at 
the  radar  receiver.  From  the  plot  it  can  be  seen  that  the  angular  error  sometimes 
exceeds  the  wingspan  of  the  aircraft,  but  usually  is  smaller  than  the  wingspan.  The 
angular  error  does  rapidly  change,  as  expected  [14].  In  order  to  see  more  detail  of 
the  scattered  field  and  the  angular  error,  “blow-ups”  of  the  two  plots  are  shown  in 
Figure  4.4.  In  this  figure,  the  quantities  are  shown  for  time  t  =  0  to  2  seconds. 

Figure  4.5  shows  the  Doppler  frequency,  fj  as  defined  in  Equation  (3.24),  and  the 
received  power  at  the  radar,  Pr  as  defined  in  Equation  (3.1),  both  versus  time.  As 
with  the  scattered  field,  “blow-up”  plots  are  shown  in  Figure  4.6.  A  radar  receiver 
must  be  able  to  track  the  changing  frequency  as  a  function  of  time.  Recall,  the  CW 
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frequency  received  by  a  radar,  fr,  is  equal  to  the  Doppler  frequency,  plus  the 
transmitter  frequency,  /.  Thus,  /,.  =  /  +  f^.  The  Doppler  bandwidth  is  a  term 
which  refers  to  the  range  of  frequencies  over  which  /,  ranges.  Since  /  does  not 
change  with  time,  the  Doppler  bandwidth  is  determined  by  the  Doppler  shift,  f^. 
This  simulation  shows  that  the  Doppler  bandwidth  is  relatively  large,  with  a  range 
of  10  kHz.  Examining  the  “blow-up”  plot  in  Figure  4.6,  the  Doppler  bandwidth  is  in 
the  range  of  300  to  3000  Hz  for  times  durations  up  to  0.1  seconds.  Obviously,  fj,  is  a 
function  of  time,  and  over  a  long  time,  fj,  may  change  by  many  kilohertz,  but  for  short 
time  durations,  under  0.1  seconds,  the  change  in  does  not  exceed  2000  Hz.  At  the 
times  when  fd  has  large  jumps,  this  usually  can  be  attributed  to  low  scattered  field 
values.  This  is  verified  by  comparing  Figures  4.3  and  4.5.  The  Doppler  bandwidths 
observed  here  seen  to  be  rather  larger  for  the  scenario,  which  might  suggest  more 
experimentation  with  the  rigid-body  and  individual  scattering  center  motion  model 
parameters  (magnitude,  cut-off  frequency,  and  correlation  matrices). 

Obviously,  only  one  simulation  example  has  been  show  here.  As  the  dynamic 
situation  changes,  different  phenomena  may  arise.  Also,  turning  off  or  changing 
the  various  aircraft  vibration  parameters  will  significantly  effect  the  scattered  field, 
angular  error,  and  Doppler  frequency. 
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-axis,  meters 


Path  of  Target 


Figure  4.1:  Path  of  C-29  aircraft  for  sample  simulation  in  {x,y,z)  space.  The  C-29 
is  flying  straight  and  level,  with  vibrations  as  described  in  the  text,  along  the  line 
shown  in  the  plot.  The  radar  is  located  at  the  origin,  denoted  by  a 


Azimuth  Angle  (AZ)  vs.  Time 


Radar-Target  Range  vs.  Time 


Elevation  Angle  (EL)  vs.  Time 


Data  Set  Used  vs.  Time 


Time,  seconds 


Figure  4.2:  Azimuth  and  elevation  angles  for  sample  simulation,  as  defined  in  Fig¬ 
ure  2.2,  versus  time,  relative  to  the  aircraft  pointing  vector  Pt.  Also,  radar-target 
range  versus  time  and  C-29  Model  Number  used  versus  time. 
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Angle  Error,  Degrees 


Figure  4.3:  Scattered  Electric  Field  and  Angular  Error  versus  time  for  sample  simula¬ 
tion.  The  dashed  lines  in  the  angular  error  plot  denote  the  relative  angular  wingspan 
of  the  aircraft  as  seen  by  the  radar  receiver. 
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Figure  4.4:  Blow-up  of  Scattered  Electric  Field  and  Angular  Error  versus  time  for  a 
portion  of  the  sample  simulation.  The  dashed  lines  in  the  angular  error  plot  denote 
the  relative  angular  wingspan  of  the  aircraft  as  seen  by  the  radar  receiver. 
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Time,  seconds 


Figure  4.6:  Blow-up  of  Doppler  frequency  and  Received  Power  versus  time  for  a 
portion  of  the  sample  simulation. 


5.  MATLAB  Code  Descriptions 


In  this  section,  all  of  the  MATLAB  codes  used  in  the  simulation  package  are 
briefly  described.  Some  of  the  codes  are  listed  in  Appendix  B.  Table  5.1  lists  all  of 
the  MATLAB  codes  needed  to  run  the  simulations.  In  Table  5.1,  the  code’s  name, 
a  brief  description  of  the  code,  and  the  type  of  MATLAB  code  is  specifled  for  each 
code.  The  two  type  of  MATLAB  codes  are  script  and  function.  A  script  code  can 
be  run  from  the  MATLAB  command  line  or  as  a  background  job.  A  function  code 
needs  to  be  properly  called  by  a  script  function  or  from  the  command  line,  and  is 
analogous  to  a  subroutine.  For  more  information  concerning  MATLAB,  consult  [18]. 
The  simulations  were  performed  using  MATLAB  version  4.2c,  dated  December  31, 
1994,  SUN  UNIX.  MATLAB  also  has  versions  which  run  on  IBM  PCs  and  APPLE 
machines. 

As  stated  before,  to  run  a  dynamic  simulation,  the  only  code  that  needs  to  be  ex¬ 
ecuted  in  MATLAB  is  profile.m.  This  script  file  calls  all  of  the  other  required  codes. 
Also,  all  options  are  definable  in  profile.m,  which  passes  the  required  information  to 

Table  5.1:  Listing  of  MATLAB  codes. 


Code  Name 

Description 

Type 

profile.m 

runs  entire  dynamic  simulation 

script 

scattered-field,  m 

calcs  various  outputs 

function 

rigid_motion.m 

calcs  rigid  body  motion 

function 

indiv_mot.m 

calcs  individual  motion 

function 

model-image. m 

plots  a  model  ISAR  image 

script 

C29model.m 

holds  all  C29  models 

function 

makemodel.m 

makes  C29  models 

script 

prony2d.m 

runs  2-D  TLS  Prony  Technique 

function 

C29image.m 

makes  C-29  ISAR  image 

script 
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the  other  MATLAB  codes.  The  script  file  profile.m  is  commented  and  the  variable 
names  have  been  chosen  to  minimize  confusion  and  to  correlate  with  the  equations 
in  this  report  as  closely  as  possible. 

The  CPU  time  of  the  simulation  performed  in  Section  4  is  on  the  order  of  two 
minutes  being  run  on  a  SUN  SPARC  2  with  32  MB  of  RAM.  Little  effort  was  given 
to  speed-up  the  simulations,  but  since  many  loops  are  in  the  MATLAB  code,  if  the 
simulations  were  rewritten  in  a  faster  language,  such  as  “C”,  then  the  required  CPU 
time  should  decrease  significantly.  MATLAB  is  an  “interpreter”  language,  and  is  not 
efficient  at  performing  loops. 

There  should  be  no  problem  converting  the  MATLAB  codes  to  another  language 
if  real-time  operation  is  required.  Note  that  it  will  only  be  necessary  to  convert  the 
following  codes  from  MATLAB  to  the  new  language;  profile.m,  scatteredJield.m, 
rigid_motion.m,  indiv_mot.m,  C29model.m,  and  modelJmage.m.  The  2-D 
TLS-Prony  Technique  codes  do  not  need  to  be  converted  as  they  are  not  used  in 
the  simulations,  they  only  create  the  models.  The  only  built-in  MATLAB  functions 
which  are  used  in  the  simulations  which  are  not  “simple”  functions  easily  performed 
in  any  language,  are  the  filter  design  functions,  firl.m,  used  in  rigid_motion.m  and 
indiv_motion.m.  The  filter  designs  can  be  performed  in  MATLAB,  and  the  results 
(filter  coefficients)  can  be  hard-coded  into  any  new  code. 


6.  Conclusions  and  Recommendations 

The  goal  of  this  simulation  package  is  to  provide  realistic  scattered  electric  field, 
Doppler  shift,  and  angular  error  values  for  a  dynamic  radar-target  engagement.  These 
values  can  be  used  in  a  hardware-in-the-loop  simulation  of  a  radar-target  or  missile- 
target  encounter  to  help  evaluate  the  performance  of  radar  receivers  or  missile  seekers. 
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It  must  be  emphasized  that  a  mono-static  radar  or  fully-active  missile  seeker  has  been 
assumed  as  the  radar  in  the  scattering  models.  However,  if  bi-static  scattering  data 
were  available,  then  the  simulations  can  be  corrected  for  this  situation. 

The  hardware-in-the-loop  simulation  would  have  to  simulate  the  scattered  electric 
field,  the  Doppler  shift,  and  the  angular  error  values  predicted  by  the  simulations. 
Details  of  how  this  is  accomplished  are  given  in  [16].  With  the  resources  available  to 
DSTO,  this  should  easily  be  accomplished. 

One  phenomena  which  is  not  taken  into  account  in  the  simulations  presented  here 
are  Jet  Engine  Modulation  (JEM)  effects  on  the  scattered  field.  JEM  effects  occur 
due  to  the  turbine  blades  in  the  target  aircraft’s  jet  (or  turbo-prop)  engines  rotating 
at  high  rates  of  speed.  JEM  effects  will  change  the  scattered  field,  angular  error,  and 
Doppler  shift  values.  For  nose-on  scenarios  (aircraft  heading  directly  toward  radar), 
the  JEM  effects  will  be  most  significant.  For  other  radar-aircraft  aspect  angles,  the 
JEM  effects  will  be  less  significant.  The  JEM  effects  tend  to  modulate  the  scattered 
field,  and  possibly  these  effects  can  be  modeled  if  the  simulations  are  upgraded.  Even 
though  the  simulations  do  not  take  JEM  effects  into  account,  the  results  are  still  valid 
and  the  JEM  effects  are  secondary  in  most  scenarios  to  the  more  important  effects  of 
aircraft  motion  (rigid  body  and  individual). 

Another  limitation  of  the  simulations  is  in  the  individual  motion  algorithm  in 
Section  2.5.2.  This  algorithm  moves  the  scattering  centers  around  randomly  about 
their  unperturbed  positions  without  any  regard  to  the  scattering  centers  location  on 
the  aircraft.  As  a  possible  upgrade  to  the  simulations,  some  weighting  factor  as  a 
function  of  scattering  center  location  could  be  added  to  the  motion  values  allowing 
scattering  centers  further  away  from  the  phase  center  of  the  target  to  move  more  than 
scattering  centers  which  are  near  the  phase  center  of  the  target.  This  is  meant  to 
simulate  scattering  centers  on  wing-tips  moving  more  than  scattering  centers  on  the 
fuselage. 
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In  the  simulations  it  is  assumed  that  the  target  (C-29  aircraft)  is  flying  and  vi¬ 
brating,  and  the  radar  is  stationary  at  position  Ro.  If  the  radar  (or  missile  seeker)  is 
also  traveling  at  a  velocity  and  vibrating,  the  simulations  can  be  upgraded  to  account 
for  this  by  varying  Ro  and  also  accounting  for  bore-sight  angle  errors. 

Taking  the  above  concerns  into  account,  the  simulations  presented  contain  a  very 
accurate  model,  the  2-D  Damped  Exponential  Model,  whose  parameters  were  derived 
directly  from  ISAR  data,  for  the  scattering  from  a  radar  target.  The  signatures  shown 
in  Section  4,  when  compared  with  actual  dynamic  signatures  [19],  show  that  there  is 
an  overall  correspondence  between  the  results  of  the  simulations  and  what  is  observed 
in  reality.  Thus,  it  is  concluded  that  this  simulation  package  would  present  actual 
radar  or  missile  receiver  hardware  with  signatures  that  are  realistic  for  a  dynamic 
radar-target  scenario. 
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8.  Appendices 


Appendix  A.  Radar  Imaging  Using  Wide- Bandwidth  Data 

This  appendix  briefly  describes  how  an  Inverse  Synthetic  Aperture  Radar  (ISAR) 
Image  is  formed  when  multiple  frequency  multiple  angle  ISAR  data  is  available.  The 
analysis  of  2-D  multiple  frequency  and  multiple  angle  radar  data  has  mainly  utilized 
the  2-D  Fourier  Transform  [20].  These  techniques  are  also  applied  in  tomography 
and  have  applications  in  electron  microscopy,  x-rays,  radio  astronomy,  geology  and 
medical  imaging. 

The  data  analyzed  in  this  effort  is  ISAR  data,  however.  Synthetic  Aperture  Radar 
(SAR)  data  could  also  be  analyzed.  The  difference  between  a  SAR  scenario  and  an 
ISAR  scenario  is  that,  in  a  SAR  scenario  the  target  is  stationary  and  the  radar  is 
moving  while  in  an  ISAR  scenario,  the  radar  is  stationary  and  the  target  is  moving.  In 
either  case,  when  a  measurement  is  taken  it  is  assumed  that  both  radar  and  target  are 
stationary.  Any  motion  present  during  a  measurement  must  be  “subtracted”  out  [8]. 
Thus,  it  is  assumed  in  this  Appendix  that  the  measurements  are  taken  are  done  so 
for  a  static  (no-motion)  situation. 

First,  the  relationship  between  the  frequency  domain  and  the  time  domain  will 
be  established.  The  2-D  image  (SAR  or  ISAR)  of  a  radar  target  is  sometimes  called 
“the  reflectivity  density  function”  [8]  and  is  denoted  h  {x,  y)  in  this  appendix. 

The  radar-target  geometry  for  2-D  imaging  purposes  is  shown  in  Figure  8.1.  The 
ttu-coordinates  are  fixed  to  the  radar  while  the  ajy- coordinates  are  fixed  to  the  target. 
The  range  between  the  maximum  extent  of  the  target  and  the  radar  is  assumed  to 
be  much  greater  than  the  maximum  extent  of  the  target.  Also,  the  target  is  assumed 
to  be  in  the  far  field  of  the  radar.  When  measurements  are  taken  by  the  radar,  it 
is  assumed  that  there  is  no  relative  velocity  between  the  radar  and  the  target.  As  a 
start,  it  is  assumed  that  the  radar  is  radiating  a  single  frequency,  /.  For  a  fixed  ip 
{ip  in  this  appendix  and  9  used  in  the  main  body  of  the  report  are  equivalent  angular 
variables)  and  a  particular  range  v,  the  signal  received  by  the  radar  is 

/OO 

h^{u,v)du  (8.1) 

-OO 
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V 


stationary  radar 


Figure  8.1:  Radar-target  geometry  for  2-D  imaging. 

where 

h^{u,v)  =  h{x,y)  (8.2) 

and 

X  =  Mcos  (■0)  —  t)  sin  ('0) 
y  =  «sin  (■0) -h  r;cos('0). 

The  total  signal  received  at  a  given  ip,  H  {ip),  is  an  integral  of  all  signals  along  the 
projection  of  h  {x,  y)  onto  the  n-axis  modified  by  a  phase  factor  which  accounts  for 
the  round-trip  phase  delay  along  the  t>-axis.  Thus, 

/oo  -i4»D 

p{v,ip)e  A  dv.  (8.3) 

-oo  ' 
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Utilizing  Equation  (8.1),  Equation  (8.3)  becomes 


r  r  h^{u,v)e^ 

J~oo  J—OO 


du  dv. 


Realizing  v  =  ycos('0)  —  a;  sin  (ip),  using  Equation  (8.2)  and  defining 


/.  = 

fv  = 


2  sin  (V') 

A 

— 2cos(V’) 
A 


the  received  signal  at  the  radar  becomes 


fOO  poo 

H{f.Jy)=  /  h{x,  y)  dx  dy. 

J  —  OO  J  ~oo 


(8.4) 


(8.5) 

(8.6) 


(8.7) 


This  is  the  standard  2-D  Fourier  Transform.  All  of  the  useful  properties  which  exist 
for  the  Fourier  Transform  exist  for  the  2-D  Fourier  Transform  pair 


h[x,y)  •< — >  H  {fxify)  •  (^•®) 

From  Fourier  Theory  [21], 

/OO  poo 

/  df,  dfy.  (8.9) 

-OO  J  —  oo 

This  relationship  can  be  written  in  polar  form  by  defining  /  =  +  /|-  Thus 

poo  p2'K 

h{x,y)  =  J^ 

relates  the  reflectivity  density  function  to  the  frequency  domain  data.  Equations  (8.1) 
through  (8.10)  form  the  Projection-Slice  Theorem  [20]. 

In  practice,  we  will  have  samples  of  H  (/,  ■0),  as  in  Figure  8.3,  and  wish  to  approx¬ 
imate  h  {x,  y)  from  these  samples.  The  samples  of  H  (/,  tp)  are  usually  the  scattering 
coefficients  of  the  target  as  a  function  of  frequency  and  angle  {S{f,6)).  Usually  one 
angular  dimension  is  fixed  and  the  other  is  varied  to  yield  a  2-D  “slice”  of  the  target. 
However,  any  angular  slice  of  the  target  can  be  used  to  generate  an  image. 

There  are  several  methods  which  are  used  to  approximate  h  {x,  y)  from  the  samples 
of  H  (/,  Ip).  If  the  samples  of  H  (/,  'tp)  are  on  a  polar  grid,  as  in  Figure  8.3,  one  method 
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Figure  8.2:  2-D  frequency  domain  data  set  for  full  angular  coverage. 

is  to  interpolate  the  data  to  a  rectangular  grid  and  apply  the  2-D  Inverse  Fast  Fourier 
Transform  (2-D  IFFT).  Also,  if  the  polar  grid  is  “almost”  rectangular,  no  interpolation 
is  performed  and  it  is  assumed  the  polar  grid  data  lies  on  a  rectangular  grid.  The 
2-D  IFFT  is  directly  applied  to  the  polar  grid  data.  This  is  the  case  for  all  the 
images  appearing  in  this  report  (Figures  2.5,  2.6,  2.7,  and  2.8).  Another  method 
is  to  numerically  approximate  the  double  integral  in  Equation  (8.10).  Yet  another 
option  is  to  use  the  convolution-backprojection  method  [20].  One  type  of  data  set 
in  the  frequency  domain  is  shown  in  Figure  8.2,  which  corresponds  to  full  angular 
coverage.  Figure  8.3  shows  the  case  of  the  data  available  for  limited  angular  coverage. 
In  practice,  data  is  usually  only  available  over  a  limited  angular  extent  (usually  the 
maximum  angular  extent  is  between  3°  and  20°  for  spotlight  mode  SAR  [8]).  For  the 
data  associated  with  the  C-29  scale  model,  the  angular  extent  is  10°  and  the  angular 
swath  is  in  the  AZ  direction. 

For  the  SAR  images  in  Figures  2.5,  2.6,  2.7,  and  2.8,  the  radar  image  was  formed 
using  the  2-D  IFFT  with  no  interpolation.  The  de-focusing  effects  of  the  polar  grid 
data  were  not  degrading  enough  to  require  focusing  (interpolation  of  2-D  frequency 
domain  data  from  polar  grid  to  rectangular  grid).  Also,  the  2-D  Damped  Exponential 
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Appendix  B.  MATLAB  Code  Listings 


The  following  codes  are  listed  in  this  Appendix:  profile.m,  scatteredJield.m, 
rigid_motion.m,  indiv_mot.m,  and  a  portion  of  C29model.ni.  Descriptions  of 
these  codes  are  given  in  Section  5. 


profile  .m 


%  MATLAB  script  file  profile.m 

%  Written  by  J.  Sacchini,  DSTO,  Adelaide,  1  May  96 
%  See  final  report  for  explanation 
% 

y,  Run  a  Typical  profile 

y. 

y  All  angles  in  degrees,  all  frequencies  in  hertz,  all  distances  in 
‘/.meters,  powers  in  watts,  gains  in  numeric  (not  dB) 

‘/,  To  turn  on  a  particular  feature,  make  the  *_on  variable  =  1,  to  turn 
‘/,  off  set  it  to  0  (actually,  anything  not  equal  to  1). 


clear  ‘/.Initially  clear  out  all  variables 


‘/  Initially  produce  SAR  image,  if  desired,  must  guess  at  AZ  and  EL  angles 
iimage=0 ; 
if  iimage==l, 
f=26.0023*10“9/3; 

AZim=-3; 

ELim=3 ; 

[px , py , a , azmid , el , f min , df , nf , da , na , choice] =C29model ( AZim , ELim , f ) ; 

"/.  [px ,  py ,  a ,  azmid ,  el ,  f min ,  df ,  nf ,  da ,  na ,  model.no]  =point .model  (AZdeg ,  ELdeg ,  f ) ; 
model. image (px , py , a , azmid , e 1 , f min , df , nf , da , na , f ) ; 
hold  on 
end 


‘/,  Input  desired  frequency,  transmitter  power  and  antenna  gain,  NOT  in  dB 
f=26.0023*10*9/3; 

Ptr=100; 

Gtr=1000; 

‘/,  Do  you  wish  warnings  of  using  models  outside  of  their  intended 
‘/,  angular  areas  and  do  you  wish  the  scattering  centers  to  be  plotted 
warnings=l ; 


59 


scplot=0; 


"/,  Input  radar  position  (fixed  for  now) 

Ro=[0  0  0]; 

%  Start  profile,  ti=initial  time,  tv=time  vector, 

%  move  target,  keep  observation  point  fixed  (for  now) 

ti=0; 

dt=0.01; 

tv= [0 : dt : 6] ; 

'/otv=0 ; 
tv=tv( : ) ; 

iterations=length(tv) ; 

In  Calculate  flight  path  positions,  without  any  aircraft  translational 

%  or  rotational  shaking 

Vt=150;  %  velocity  of  target 

D't=[-1  0  0];  %  direction  of  target 

Dt=Dt/norm(Dt) ;  */,  normalize  direction 

Rtinit=[2800  -100  -100];  '/oinitial  target  position 

Rtv_nos=(Rtinit’*ones(l,iterations)) ’+Vt*(Dt’*tv’) ' ; 

Dtv_nos=(Dt'*ones(l, iterations)) ’ ; 

%  Put  in  fixed  roll,  pitch,  and  yaw  angles  if  desired 
Roll_Angle=0;  %  Constant  Roll  Angle 
Pitch_Angle=0;  %  Angle  of  Attack 
Yaw_Angle=0;  %  Crab  Angle 

RPYv_nos=[Roll_Angle*ones (iterations,!)  Pitch_Angle* . . . 
ones (iterations,!)  Yaw_Angle*ones (iterations , 1)]  ; 

/  Translational  and  rotational  motion  shaking  (rigid  aircraft) 
randn( ' seed’ , 1)  %  start  with  same  seed  each  time  if  desired 
trans_on=l ; 

%  Cut-off  freq  of  translational  shaking  in  Hz 
mag_trans=0 . 1 ;  %  Standard  Deviation  of  Magnitude  of  shaking  in  meters 

rot_on=l ; 

fc_rot=2;  %  Cut-off  freq  of  rotational  shaking  in  Hz 
niag_rot_deg=2 ;  j.  Standard  Deviation  of  Angle  shaking  in  degrees 
[Rtv_s,Dtv_s,RPYv_s]=rigid_motion(Rtv_nos,Dtv_nos,RPYv_nos,tv, . . . 
dt , trans.on , f  c.trans ,mag_trans , rot.on , f c_rot ,mag_rot_deg) ; 

Rtv=Rtv_s;  ^position 
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Dtv-Dtv_s;  '/odirection  of  flight,  not  implemented 
RPYv=RPYv_s;  %roll, pitch, yaw  in  degrees 

lnltUt/tUtlt  Insert  other  rigid  body  motion  models  here  -  Dr.  Mark  Petrusma 
y,  make  Rtv  be  the  x,y,z  perturbations  in  meters,  Rtv  has 
%  iterations  rows  and  3  columns 

y.  Make  RPYv  be  the  roll,  pitch,  yaw  perturbations  in  degrees,  RPYv  has 
y,  iterations  rows  and  3  columns 

y.  Make  sure  if  you  use  alternate  models,  to  turn  off  trans  and  rot 
%  motion  above 

y.  Individual  Scattering  Center  Motion  (Flexing  Aircraft) 
indiv_mot_on=l ; 
f c_indiv=10; 
mag_indiv_mot=0 . 02 ; 

nkeep=10;  %  must  specify  how  many  scattering  centers  to  move  independently 
IMvx=zeros (iterations, nkeep) ; 

IMvy=zeros(iterations,nkeep) ; 
for  k=l: nkeep, 

IMvx ( : , k) =indiv_mot (tv , dt , indiv.mot _on , f c_ indiv ,mag_indiv_mot ) ; 

IMvy ( : , k)=indiv_mot (tv , dt , indiv_mot_on , f c_indiv ,mag_indiv_mot ) ; 
end 

/  Now  run  profile  for  all  times  in  tv  and  calculate  all  quantities 

y. 

for  ti=l : iterations , 
t=tv(ti) ; 

Rt=Rtv(ti, : ) ; 

Rt_nos=Rtv_nos (ti , : ) ; 

Dt=Dtv(ti, :) ; 

RPY=RPYv(ti,:); 

if  ti==l,  Rt_delay=Rtv(l, :) ;  vt_old=0;  Rsc_old=0;  end 
if  ti~=l,  Rt_delay=Rtv(ti-l, :) ;  end 
y,  Dt=-(Rt-Ro);  %  to  point  directly  at  obs  point 
aRto (ti)=norm(Rt~Ro) ;  %  Range  between  target  and  observation  point 
[Es(ti) ,Es_complex(ti) ,D0A(ti) ,Pr(ti) ,fd(ti) ,AZdeg(ti) ,ELdeg(ti) , . . . 
azmid(ti) ,el(ti) ,model_no(ti) ,vt_old,Rsc_old,pha_old,extent_deg(ti) . . . 
,dvec(ti, :) ,d(ti, :)]=scattered_f ield(t,f ,Ro,Rt,Rt_nos,vt_old, . . . 

Rsc_old ,pha_old, dt ,Dt ,RPY , Vt ,Ptr ,Gtr , warnings , scplot , ti , IMvx , . . . 

IMvy , indiv_mot_on) ; 
end 
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if  iterations>l,  fd(l)=fd(2);  end 

%  Look  for  places  where  we  changed  model  numbers,  and  throw  out 
Xdoppler  because  it  is  incorrect  during  model  cheingeovers . 
model_diff=model_no-[model_no(l)  model_no(l riterations-l)] ; 
ich=f ind(model_diff ~=0) ; 
fd(ich)=fd(ich-l) ; 

y,  Plot  output  data,  if  desired 
iplot=l; 
if  iplot==l, 
figure 

subplot (2, 1,1) 
plot(tv,fd) 

xlabeK ’Time,  seconds’) 

ylabeK ’Doppler  Frequency,  Hz’) 

title ( ’Doppler  Frequency  vs.  Time’) 

subplot (2, 1,2) 

plot (tv , 10*logl0 (Pr) ) 

xlabeK ’Time ,  seconds’) 

ylabeK ’Received  Power  at  Radar,  dB’) 

title ( ’Received  Power  vs.  Time’) 

figure 

subplot (2, 1,1) 

plot (tv , 20*logl0 (abs (Es ) ) ) 

xlabeK ’Time,  seconds’) 

ylabeK ’Magnitude  of  E-Field,  dB’) 

title ( ’Scattered  Electric  Field  at  Receiver  vs.  Time’) 

subplot (2, 1,2) 

plot(tv,D0A) 

xlabel( ’Time ,  seconds’) 

ylabeK ’Angle  Error,  Degrees’) 

title (’Angle  Error  vs.  Time’) 

hold  on 

plot (tv , extent _deg/2 , ’  —  ’ ) 
plot (tv , -extent_deg/2 , ’ — ’ ) 
hold  off 
figure 

subplot (2,2,1) 
plot(tv,AZdeg) 

%  hold  on 

%  plot (t V , azmid ,’-.’) 
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“/o  hold  off 

xlabeK 'Time,  seconds’) 
ylabel ( ’ Azimuth ,  degrees ’ ) 
title( ’Azimuth  Angle  vs.  Time’) 
subplot (2, 2, 2) 
plot(tv,ELdeg) 

‘/o  hold  on 
%  plot (tv, el, ’-. ’) 

“/.  hold  off 

xlabelC’Time,  seconds’) 
ylabel ( ’Elevation  Angle,  Degrees’) 
title ( ’Elevation  Angle  vs.  Time’) 
y,  figure 
subplot (2, 2, 3) 
plot(tv,aRto) 
xlabeK ’Time ,  seconds’) 
ylabel ( ’Radar-Target  Range,  meters’) 
title ( ’Radar-Target  Reinge  vs.  Time’) 
subplot (2,2,4) 
plot (tv ,model_no) 
xlabeK ’Time ,  seconds’) 
ylabel (’C29  Model  Number’) 
title (’Data  Set  Used  vs.  Time’) 
figure 

plot3(Rtv( : ,1) ,Rtv( : ,3) ,Rtv( : ,2)) 

view([l  1  1]); 

grid 

xlabeK ’x-axis,  meters’) 
ylabeK ’z-axis ,  meters’) 
zlabeK ’y-axis,  meters’) 
title (’Path  of  Target’) 
hold  on 

plot3(Ro(l),Ro(3),Ro(2), ’*’) 
hold  off 
end 
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%  %*/.'/.'/.ny.'/.‘my,y.  scattered.f  ield  .m 


function  [Es,Es_complex,DOA,Pr,fd,AZdeg,ELdeg,azmid,el,model_no, . . . 
vt_old,Rsc_old,pha_old,extent_deg,dvec,d]=. . . 

scattered.f ield (t , f , Ro , Rt , Rt _nos , vt_old , Rsc_old , pha.old , dt , Dt , . . . 

RPY , Vt , Ptr , Gtr , warnings , scplot , t i , IMvx , IMvy , indi v_mot_on) ; 

y,  This  function  calculates  the  scattered  field,  Es,  in  volts/meter, 

%  and  the  angular  error,  in  degrees,  for  the  radar 
y,  scattering  from  a  target,  as  described  below.  Also  calculated  are 
y  the  doppler  frequency,  fd,  and  the  average  power  received  at  the 
y,  radar,  Pr.  The  input  parameters  are  the  CW  frequency  f, 
y,  the  observation  point  Ro=[xo,yo,zo] ,  the  target  location,  centered  at 
y,  Rt=  [xt , yt , zt]  and  traveling  in  the  direction  of  the  direction  vector 
%  of  Dt=[vx,vy,vz] .  The  target  can  be  traveling  with  a  velocity  of  Vt 
%  in  the  Dt  direction.  Note  that  Vt  should  be  a  scalar, 
y  It  is  assumed  that  the  time  dependence,  e^Cjut)  is  suppressed  for 
y  the  calculation  of  Es_complex,  but  is  included  for  Es  calculation, 
y  Thus,  the  current  time,  t,  should  be  input  in  seconds, 
y  The  target  is  described  by  the  damped  exponentials  (px,py,a) 
y,  which  represent  the  2-D  spectral  domain  signature  of  the  target.  I 

y  have  assumed  that  the  target  is  parallel  to  the  xz-plane  and  pointing 
y  in  the  direction  described  by  Pt.  Pt  is  defined  as  the  waterline  of 
y  the  aircraft.  Pt  and  Dt  may  differ  if  roll,  pitch,  or  yaw  are 
y  present.  The  power  of  the  transmitter,  Ptr,  and  the  antenna 
y  transmitter  gain,  Gtr,  must  be  supplied.  All  of  the  info  concerning 
y  the  radar  (or  scattering)  data  the  model  is  based  on  is  given  in 
y  C29model.m,  i.e.  fmin,df ,nf ,amin,da,na,  where 
y  fmin=lower  frequency  of  radar  data,  in  Hz 
y  df=frequency  step  of  radar  data,  in  Hz 
y  nf=number  of  frequencies  of  radar  data,  integer 
y  azmid=middle  azimuth  angle  for  input  data  between  centerline  of 
y  aircraft  and  center  of  radar  data,  in  degrees 

y  da=angular  step,  in  degrees 

y  na=number  of  angular  data  points  of  radar  data,  integer 
y  el=elevation  of  data  set,  in  degrees 

y 

y  Author:  Joe  Sacchini,  jsacchin@afit.af.mil,  1  May  96,  Australia 
y  Please  see  accompanying  final  report  for  more  detail 

y  Calculate  AZ  and  EL  angle ,  as  seen  by  the  target ,  from  target  to 


64 


%  observation  point 

Rot=Rt-Ro;  7o  Observation  point  to  Target  distance 

y,  AZ  and  EL  calculations,  center  of  target,  traveling  in  Dt  dir,  to 
"/oobservation  point,  target  is  traveling  in  Dt  dir,  but  can  have  ROLL, 
'/.PITCH,  or  YAW,  RPY  correction 

'/.  First  AZ  calculation  assuming  no  roll,  pitch,  or  yaw 
Rotaz=[Rot(l)  0  Rot(3)]; 

Dtaz=[Dt(l)  0  Dt(3)] ; 
if  norm(Rotaz) ~=0, 

AZrad=acos (sum( (-Rotaz) . *Dtaz) / (norm(Rotaz) *norm(Dtaz) ) ) ; 
end 

if  norm(Rotaz)==0,  AZrad=0;  end 
AZsign=cross (-Rotaz, Dtaz) ; 

AZrad=AZrad*sign(AZsign(2) ) ; 

'/.  Next,  EL  calculation  assuming  no  roll,  pitch,  or  yaw 
Rotel=[Rot(l)  Rot(2)  0]; 

Dtel=[Dt(l)  Dt(2)  0]; 
if  norm(Rotel)~=0, 

ELrad=acos (sum(Rotel . * (-Dtel) ) / (norm(Rotel) *norm(Dtel) ) ) ; 
end 

if  norm(Rotel)==0,  ELrad=0;  end 
ELsign=cross(-Rotel,Dtel) ; 

ELrad=ELrad*sign(ELsign(3) ) ; 

'/.  Now  correct  AZ  and  EL  for  ROLL,  PITCH,  and  YAW 
R0LL=RPY(1) ♦pi/180;  PITCH=RPY(2) ♦pi/180;  YAW=RPY(3)^pi/180; 
AZradn=AZrad+YAW ; 

ELradn=ELrad+PITCH ; 

AZradn=atan2(-sin(ELradn)^sin(R0LL)+sin(AZradn)^cos(ELradn)^cos(R0LL) , . . 

cos(AZradn)^cos(ELradn)) ; 

ELradn=asin(sin(ELradn)^cos (ROLL)+sin(AZradn) ♦cos (ELradn) ♦sin (ROLL) ) ; 

'/,  Convert  to  degrees 
AZdeg=AZradn^l80/pi ; 

ELdeg=ELradn^l80/pi ; 

'/,  Determine  physical  extent  of  aircraft  for  given  AZ  and  EL  angles, 
'/.assuming  an  C29  aircraft  fits  in  a  box  of  dimensions  SD=(SX,SY,SZ) 
'/.meters  and  it  is  Rot  meters  away  from  the  observation  point 


65 


SD=[15.6  5.37  15.6]; 

extent=(SD(3)*cos (AZradn)+SD(l) *sin(AZradn) * (abs (AZradn) >pi/4) ) ; 
extent_deg=(extent/nonn(Rot))*(180/pi) ; 

y,  Determine  which  C29  model  to  use 

[px , py , a , azmid , el , imin , df , nf , da , na , model.no] =C29model ( AZdeg , ELdeg , f ) ; 

% [px , py , a , azmid , el , f min , df , nf , da , na , model_no] =po int _model (AZdeg , ELdeg , f ) ; 

7.  Determine  how  to  use  model,  what  m  and  n  are 

no_scatctrs=length(px) ; 

c=2.9979*10~8; 

azmid_rad=azmid*pi/180 ; 

nAZrad=AZradn-azmid_rad ; 

nAZdeg=nAZrad*180/pi ; 

darad=da*pi/180 ; 

ko=2*pi*f/c; 

Dt=Dt/norm(Dt) ; 

Ruf=c/ (2*df ) ; 
dfa=f*darad; 

Rua=c/(2*df a) ; 
xstep=Rua/na; 
ystep=Ruf /nf ; 

Dt=Dt/norm(Dt) ; 

ELwarning=5/2 ; 
if  warnings==l, 
if  abs(el-ELdeg)  >  ELwarning, 
disp( ’Warning,  Elevation  off  by’) 
disp(el-ELdeg) 
end 
end 

%  Calculate  offset  from  center  of  grid  in  2D  frequency  space,  and 

7,  check  for  validity  of  model,  m=freq  no,  n=angle  no 

m=(f-fmin)/df ;  %  Frequency  step 

nfmin=-nf/2; 

nfmax=nf*1.5; 

if  m<nfmin,  m=nfmin;  flagmlow=l;  end 
if  m>nfmax,  m=nfmax;  flagmhigh=l;  end 
n=na/2+nAZrad/darad;  %  Angular  step 
namin=-na/2 ; 
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nainax=na*l .  5 ; 

if  iKnamin,  n=namin;  flagnlow=l;  end 
if  n>namax,  n=namax:  flagnhigh=l;  end 
if  warnings==l, 

if  f laginlow==l ,  disp( ’Warning,  frequency  too  low  for  valid  model’);  end 
if  f lagmhigh==l ,  disp( ’Warning,  frequency  too  high  for  valid  model’);  end 
if  flagnlow==l,  disp( ’Warning,  azimuth  too  low  for  valid  model’);  end 
if  f lagnhigh==l ,  disp( ’Warning,  azimuth  too  high  for  valid  model’);  end 
end 
y.[m,n] 

%  Loop  through  each  scattering  center 
for  k=l :no_scatctrs , 

%  Randomly  displace  the  scatting  center,  if  desired,  independent  of 
"/,  the  other  scattering  centers  (individual  motion  or  non-rigid  body  motion) 
indiv_corr=[0.8  0.2;  0.2  0.8]; 
if  indiv_mot_on==l , 
x_mot=IMvx(ti ,k) ; 
y_mot=IMvy (ti , k) ; 

corr_mot=indiv_corr* [x_mot ; y_mot] ; 
x_mot=corr_mot(l) ; 
y_mot=corr_mot (2) ; 
else 
x_mot=0 ; 
y_mot=0  ; 
end 

px (k) =abs (px (k) ) ♦exp ( j  * (angle (px (k) ) +x_mot ) ) ; 
py (k) =abs (py (k) ) ♦exp ( j ♦ (angle (py (k) ) +y_mot ) ) ; 
if  scplot==l, 

plot (angle (py(k))^Rua/(2^pi)-0.5^xstep, . . . 

-angle (px (k) ) ♦Ruf / (2^pi) +0 . S^ystep , ’ + ’ ) 

hold  on 
end 

%  Calcuate  the  range  to  the  k-th  scattering  center,  in  3-D,  for  DOA 
%  estimation,  azrot  is  the  angle  between  Dt  and  the  xy-pleine  while 
%  elrot  is  the  angle  between  Dt  and  the  xz-plane.  azrot  and  elrot 
%  have  no  dependence  on  absolute  target  or  observation  location,  thus 
%  they  are,  in  general,  different  than  AZ  and  EL 
Rscl=[-angle(px(k))^Ruf/(2^pi)  0  -angle (py (k) ) ♦Rua/ ( 2^pi)]  ; 
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eizrot=atan(Dt(3)/Dt(l)) ; 
elrot=atan(-Dt (2)/Dt(l) ) ; 

%  Correct  for  ROLL,  PITCH,  and  YAW 

R0LL=RPY(l)*pi/180;  PITCH=RPY(2)*pi/180;  YAW=RPY(3) ♦pi/180; 
azrotn=azrot+YAW ; 
elrotn=elrot+PITCH ; 

azrotn=atan2(-sin(elrotn)^sin(R0LL)+sin(azrotn)*cos(elrotn)*cos(R0LL) , . . 

cos(azrotn)*cos(elrotn)) ; 

elrotn=asin (sin (elrotn)*cos (ROLL) +sin(azrotn)*cos(elrotn)*sin(ROLL) ) ; 
y,  now  create  rotation  matrix  and  determine  distance  and  direction  from 
7,  the  k-th  scattering  center  to  the  observation  point 
rot_mata=[cos(azrotn)  0  sin(azrotn) 

0  10 

-sin(azrotn)  0  cos(azrotn)]  ;  /(Rotate  in  AZ 
rot_mate=[cos(elrotn)  -sin(elrotn)  0 
sin(elrotn)  cos(elrotn)  0 

0  0  1]  ;  y,Rotate  in  EL 

rot_mat=rot_mata^rot_mate ;  %  Total  rotation  matrix 
Rsc(k, :  )=(rot_mat+Rscl’ ) '+Rt;  yodistance  from  k-th  SC  to  origin 
Rsco(k, : )=Rsc(k, : )-Ro;  ‘/(distance  from  k-th  SC  to  obs  point 
vt(k, :)=-Rsco(k, :)/norm(Rsco(k, :)) ;  %  direction  from  k-th  SC  to  obs  pt 

‘/,  Calculate  damping  factor  inherent  in  damped  exponential  model  (may 
*/,be  turned  off  if  desired).  If  this  is  turned  off,  then  all  you  have 
‘/,is  a  true  point  scatterer  model.  Also,  the  damping  has  been  limited 
‘/,to  a  factor  of  1000  (arbitrary) . 
no_damping=0;  ‘/,  to  turn  off  damping,  make  no_damping=l 
damping=abs (px (k) ) “m+abs (py (k) ) “n ; 
if  damping  <  0.001,  damp ing=0. 001;  end 
if  damping  >  1000,  damping=1000 ;  end 
if  no_damping==l ,  damping=l;  end 

'/,  Calculate  the  scattered  field  due  to  the  k-th  scattering  center  at 
’/the  center  of  the  target 
mag (k) =abs ( a(k) ) ♦damping ; 

pha (k) =angle (a (k) ) +angle (px (k) ) ♦m+angle (py (k) ) ♦n ; 
Es_s(k)=mag(k)^exp(j^pha(k) ) ; 

end 

‘/o  Calculate  distance  between  target  center  and  observation  pt,  r 
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r=norm(Rt-Ro) ; 

ko=2*pi*l/c; 

etao=376.7; 

k  Calculate  total  scattered  field  at  center  of  target 
Es_complex=suin(Es_s)  ; 

%  Allow  for  all  power  and  range  calculations 
if  r~=0, 

Es_complex=Es_complex*sqrt(2*Ptr*Gtr*etao)/(4*pi*r*2)*exp(-j*ko*2*r) ; 

end 

'/«  Put  in  e“{jwt}  time  dependence  and  go  from  spectral  to  time 
Es=abs (Es.complex) *cos (2*pi*f *t+angle (Es.complex) ) ; 


y.  Calculate  DOA  estimate,  R  here  is  Ro-Rt,  dir  from  target  center 
%  to  obs  point 

“/.(working  in  a  target  fixed  coord-system  for  DOA  calcs) 

*/,  NOTE:  DOA  here  means  angluar  error!!!!! 
d=[0  0  0]; 

mag=mag*sqrt(2*Ptr*Gtr*etao)/(4*pi*r‘2) ; 

R=Ro-Rt ; 

R_nos=Ro-Rt_nos ; 
for  k=l :no_scatctrs , 
aa(k)=0; 

for  1=1 :no_scatctrs, 
aa(k)=aa(k)+mag(k)*mag(l) . . . 

*cos ( (pha(k) -2=i'ko*sum(vt  (k , : )  .  *R) ) - (pha(l) -2*ko*sum(vt  (1 , : )  .  *R)  )  )  ; 
end  ’ 

d=d+aa(k)*vt(k, : ) ; 
end 

d=d*  (ko/  (abs  (Es_complex)  .  '‘2) )  ; 

D0Arad=acos ( sum(R_nos . *d) / (norm(R_nos ) ♦norm(d) ) ) ; 

DOAdir=cross(d,R) ; 
if  D0Adir(3)~=0, 

D0Asign=sign(D0Adir(3)) ; 
else 

D0Asign=sign(D0Adir(2)) ; 
end 

D0A=D0Arad*18O/pi*D0Asign ; 
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"/,  Calculate  estimate  of  received  power,  assuming  the  same  receiving 
and  transmitting  antenna 
Pr=abs (Es.complex)  '‘2* (2*pi*Gtr) /  (etao* (c/f )  "2) ; 

y,  Calculate  the  doppler  frequency,  method  one 
y,dr_dt=  (Rt-Rt  .delay )  /dt ; 
y,vel=(norm(dr_dt  ,2)) ; 
y.fd=sum(2*d.*dr_dt)/(2*pi) ; 

y,  Doppler  method  two,  better  method 

wd=0; 

if  ti"=l, 

for  k=l :no_scatctrs , 
d_phi= (pha(k) -pha.old (k) ) /dt ; 

d_dir=sum((vt(k, : )-vt_old(k, : )) .*(Rsc(k, : )))/dt ; 
d_range=sum((Rsc(k, : )-Rsc_old(k, : )) ,*vt(k, : ))/dt; 
wd=wd+aa(k) * (d_phi-2*ko* (d_range+d_dir ) ) ; 
end 
end 

dvec=[d_phi  d.range  d_dir  d_phi  -ko*(d_range+d_dir)] ; 
if  ti==l,  dvec=[0  0000];  end 
wd=wd*(-l/(abs(Es_complex) . “2)) ; 
fd=wd/ (2*pi) ; 

dopp_lim=0;  */.  only  use  if  doppler  shift  spikey  -  this  is  a  fudge  factor 
if  dopp_lim==l,  y,  limit  doppler  spikes  due  to  field  drop  out 
fdmax=1.5*((2*Vt*f)/(3*10‘8)); 
fdmin=0.5*((2*Vt*f )/(3*10*8)) ; 
if  fd  >  fdmax,  fd=fdmax;  end 
if  fd  <  fdmin,  fd=fdmin;  end 
end 

vt_old=vt ; 

Rsc_old=Rsc ; 
pha_old=pha; 
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liyiymilU  rigid_mot  ion .  m  m/.m"/.'/.'/.’/.'/.'/.'/.*/. 


function  [Rtv_s,Dtv_s,RPYv_s]=rigid_motion(Rtv_nos,Dtv_nos, . . . 

RPYv_nos,tv,dt,trans_on,fc_trans,mag_trans,rot_on,fc_rot ,mag_rot_deg) ; 

y.  Do  translational  shaking 
if  treins_on~=l , 

Rtv_s=Rtv_nos ;  %  NO  shaking 

else  %  else  shake 
nts=length(tv) ; 
fs=l/dt; 

%  correlation  of  x  y  and  z  components 
corr_mat=[0.8  0.1  0.1;0.1  0.8  O.ljO.l  0.1  0.8]; 

7,  if  cut-off  freq  is  under  Nyquist 
if  fc_trans<(f s/2) , 
f c_transn=f c_trans/f s ; 
wnx=randn(nts , 1) *mag_trans ; 
wny=randn(nts , l)*mag_trans ; 
wnz=randn(nts , 1) *mag_trans ; 
b=f irl(16,f c_trcinsn/2) ; 
wnfx=filter(b,l,wnx) ; 
wnfy=filter(b,l,wny) ; 
wnfz=f ilter(b, 1 ,wnz) ; 
wnf=[wnfx  wnfy  wnfz] ; 
else 

7.  If  cut-off  freq  exceeds  Nyquist 
f sn=8*fc_trans ; 
up=round(f sn/f s) ; 
f sn=up*f s ; 

wnx=randn (nt s*up , 1 ) *mag_trans ; 
wny=randn (nt s  *up , 1 ) ♦mag_trans ; 
wnz=randn (nts*up , 1 ) *mag_trans ; 
b=firl(16, (f c_trans/fsn)/2) ; 
wnfx=filter(b,l,wnx) ; 
wnfy=filter(b,l,wny) ; 
wnf z=f ilter (b , 1 , wnz) ; 

wnf=[wnfx(l:up:length(wnfx))  wnfy(l:up:length(wnfy)) . . . 

wnfz(l;up:length(wnfz))] ; 

end 

wnf =wnf *corr_mat ;  %  Correlate  x  y  and  z  movement 
Rtv_s=Rtv_nos+wnf ;  %  Add  shake  to  positions 
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end 


y,  Now  do  ROLL,  PITCH,  and  YAW  shaking 
if  rot_on~=l, 

RPYv_s=RPYv_nos ; 
else 

mag_rot=mag_rot_deg* (pi/180) ;  %  Magnitude  of  shaking  in  radians 

nts=length(tv) ; 
f s=l/dt ; 

%  correlation  of  roll  pitch  and  yaw  components 
corr_mat=[0 .6  0.2  0.2;0.2  0.6  0.2;0.2  0.2  0.6]; 

'/,  if  cut-off  freq  is  under  Nyquist 
if  f c_rot<(f s/2) , 
f c_rotn=f c_rot/f s ; 
wnx=raiidn(nts ,  l)*mag_rot ; 
wny=randn(nts , l)*mag_rot ; 
wnz=randn(nts , l)*mag_rot ; 
b=f irl(16,f c_rotn/2) ; 
wnfx=f ilterCb, l,wnx) ; 
wnfy=filter(b,l,wny) ; 
wnfz=f ilter(b,l,wnz) ; 
wnf=[wnfx  wnfy  wnf z] ; 
else 

%  If  cut-off  freq  exceeds  Nyquist 
f sn=8*f c_rot ; 
up=round(f sn/f s) ; 
f sn=up*f s ; 

wnx=randn (nt  s  *up , 1 ) *mag_r ot ; 
wny=randn(nts*up, l)*mag_rot ; 
wnz=randn(nts*up, l)*mag_rot ; 
b=firl(16, (f c_rot/f sn)/2) ; 
wnfx=f ilter(b,l,wnx) ; 
wnfy=f ilter(b,l,wny) ; 
wnfz=filter(b,l,wnz) ; 

wnf=[wnfx(l :up:length(wnfx))  wnfy(l :up:length(wnfy)) . . . 

wnfz (1 : up: length (wnf z))] ; 

end 

wnf=wnf*corr_mat;  /,  Correlate  ROLL  PITCH  and  YAW  movements 
RPYv_s=RPYv_nos+wnf ;  %  Add  shake  to  angles 

end 
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%  Now  do  direction  of  flight  shaking,  if  desired 
Dtv_s=Dtv_nos;  %  No  change  in  direction  implemented  yet 
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indiv_mot  .m  y.r/.7.%/x/.y, 7.7.7, 


function  IMv=indiv_mot (t v , dt , indiv_mot_on , f c.indiv ,mag_indiv_mot ) ; 

nts=length(tv) ; 

IMv_no=zeros (nts ,1) ; 


7.  Do  individual  Scattering  Center  shaking 
if  indiv_mot_on~=l , 

IMv=IMv_no;  7.  NO  shciking 
else  7.  else  shcike 
f s=l/dt ; 

7,  if  cut-off  freq  is  under  Nyquist 
if  f c_indiv<(f s/2) , 
f c_indivn=f c_indiv/f s ; 
wn=randn(nts , l)*mag_indiv_mot ; 
b=f irl(16,f c_indivn/2) ; 
wnf=filter(b,l,wn) ; 
else 

7,  If  cut-off  freq  exceeds  Nyquist 
f sn=8*f c.indivs ; 
up=round(f sn/f s) ; 
fsn=up*fs; 

wn=randn(nts*up , 1) *mag_indiv_mot ; 
b=firl(16, (f c_indiv/f sn)/2) ; 
wnf =f ilter (b , 1 , wn) ; 
wnf=wnf (1 : up: length (wnf)) ; 
end 

IMv=wnf;  7.  Add  shake  to  positions 
end 
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I  rmmm  portion  of  c29modei.m  v.y.mrmy.m 


function  [px,py,a,azmid,el,fmin,df ,nf ,da,na,model_no]=. . . 

C29model ( AZdeg , ELdeg , f , display) ; 


'/o  This  function  determines  which  C29-based  2-D  damped  exponential 
%  model  is  best  to  use  for  the  desired  AZ  and  EL  angles .  The  AZ  and  EL 
"/,  angles  are  defined  relative  to  the  target,  the  C29.  The  outputs  are 
%  the  target  model  [px,py,a] ,  the  middle  azimuth  of  the  data,  azmid 
*/,  (degrees),  the  elevation  angle  of  the  data,  el  (degrees),  the  min 
‘/,  freq,  fmin  (Hz),  the  frequency  step,  df  (Hz),  the  number  of  frequency 
’/,  samples,  the  angular  step,  da  (degrees),  and  the  number  of  eingles. 

%  This  function  will  find  the  "closest"  data  set  to  the  specified  AZ  and 
'/,  EL  angles.  The  desired  frequency,  f,  can  also  be  factored  into  the 
y,  calculation,  but  is  not  presently. 

%  Display  images,  of  original  data  and  scattering  center  locations, 
y,  and  model  generated  SAR  image  and  scattering  center  locations,  if  you 
y,  wish  by  making  display=l. 

y, 

y.  ALL  ANGLES  ARE  IN  DEGREES,  NO  EXCEPTIONS  IN  THIS  CODE 

if  nargin<4  display=0;  end 
if  narginO  f=(26.0023*10"9)/3;  end 

y,  Azimuth  and  elevation  angles  of  model  data 

AZmodel=C  -3.76199340820312 

-0.08099365234375 

3.71899414062500 

-3.76100158691406 

-0.08099365234375 

3.71899414062500 

-3.76100158691406 

-0.08099365234375 

3.71899414062500 
20.23899841308594 
23.91900634765625 

27.71899414062500 
20.23899841308594 
23.91900634765625 

27.71899414062500 
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20.23899841308594 

23.91900634765625 

27.71899414062500 

-20.23899841308594 

-23.91900634765625 

-27.71899414062500 

-20.23899841308594 

-23.91900634765625 

-27.71899414062500 

-20.23899841308594 

-23.91900634765625 

-27.71899414062500] ; 

ELmodel= [3 . 00699996948242 

3.00699996948242 

3.00699996948242 

5.00600051879883 

5.00600051879883 

5.00600051879883 

7.00900077819824 

7.00900077819824 

7.00900077819824 

3.00699996948242 

3.00699996948242 

3.00699996948242 

5.00600051879883 

5.00600051879883 

5.00600051879883 

7.00600051879883 

7.00600051879883 

7.00600051879883 

3.00699996948242 

3.00699996948242 

3.00699996948242 

5.00600051879883 

5.00600051879883 

5.00600051879883 

7.00600051879883 

7.00600051879883 

7.00600051879883] ; 


%  Determine  differences  between  desired  angles  and  model  angles 
AZdif f =AZmodel-AZdeg*ones (size (AZmodel) ) ; 

ELdif f =ELmodel-ELdeg*ones (size (ELmodel) ) ; 


%  Define  sorting  metric  as  total  angular  error 
metric=abs ( AZdif f)+abs(ELdiff) ; 

%  Sort  data  by  metric  from  actual  AZ  EL  to  available  model  AZ  EL. 
[sorted.metric , rank] =sort (metric ) ; 

%  Choose  the  AZ  and  EL  model  with  the  smallest  metric 
model_no=rank(l) ; 

%  Begin  listing  model_no  1  through  27,  and  choose  the  appropriate  model 

if  model_no==l, 

7,  DATA  SET  1 

%  C29  data,  C29_3D_hhN, 02700,  nkeep=10,  dataf=G(l :2 : 128 , 1 :64) 

7.  AZ=  -3.76199340820312 
7.  EL=  3.00699996948242 
nkeep=10; 

fmin=26.0023*10‘9/3; 

df=20*10“6/3; 

nf=64; 

da=0 . 04 ; 

na=64 ; 

azmid  =  -3.76199340820312; 
el  =  3.00699996948242; 
energy  =  l.Oe+02  *[4.46351955889678 
3.38992633084404 
3.09244961258918 
3.03737934157418 
2.11876036758480 
1.55804105492920 
1.44295351235822 
1.32096395403193 
0.57293375351862 
0.50675258067672] ; 

px  =[1.00734824451613  -  0.120209318992381 
0.34346787678126  -  1.048924288514941 
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0.19620258114083  -  1.142809946529011 
0.93987525649699  -  0.370950798569281 
0.98559136671387  +  0.200215199001371 
0.83608629828962  -  0.585169777010671 
-0.16587689351521  -  1.026620433443001 
1.00734824451613  -  0.120209318992381 
0.93987525649699  -  0.370950798569281 
-0.12059389273158  +  0.983134820363631]; 
py  =[0.92461720846197  -  0.389069078246391 
1.00169462678913  +  0.004954742794401 
1.00169462678913  +  0.004954742794401 
1.00169462678913  +  0.004954742794401 
1.00169462678913  +  0.004954742794401 
0.92461720846197  -  0.389069078246391 
1.00169462678913  +  0.004954742794401 
1.00169462678913  +  0.004954742794401 
0.95422339145316  -  0.278144597708011 
0.95773301389047  +  0.321698724883851]; 
a  =[-0.11870475204384  +  0.131064559413671 
0.00159621059295  +  0.000912966817221 
0.00009043975091  -  0.000025404784491 
0.10690730553061  +  0.143964727536751 
-0.00016535310040  -  0.177927208093691 
0.00076855880197  +  0.081245343998841 
-0.03155655411428  -  6.010359594132321 
-0.07005013611894  +  0.072577675030991 
-0.01078083345972  -  0.097798567183411 
-0.06936266338304  +  0.074463684026151]; 

re  =[0.10466510944782  0.71533402626833  0.98584320614570]; 

end 

V/tV/tV/'/tVh  26  more  models  are  listed  In  the  C29model.m,  but  not  here  "/X/X 
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